# source:git/Singular/LIB/divisors.lib@389c6d

spielwiese
Last change on this file since 389c6d was 389c6d, checked in by Janko Boehm <boehm@…>, 9 years ago
Fixed a bug in divisors.lib
• Property mode set to `100644`
File size: 22.2 KB
Line
1///////////////////////////////////////////////////////////////////////////////
2version="version divisors.lib 4.0.0.0 Jun_2013 "; // \$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 interger
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)            evalutate 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 interger
56degreeFormalDivisor(formaldivisor)              degree of a formal divisor
57
58makePDivisor(List)                              make a formal polyhedral sum of divisors
59evaluatePDivisor(pdivisor,intvec)               evaluate a polyhedral divisor
60                                                to an integer formal divisor
61pdivisorplus(pdivisor,pdivisor)                 add two polyhedral divisors
62
63
64";
65
66
67////////////////////////////////////////////////////////////////////////////////
68
69LIB "primdec.lib";
70
71proc mod_init()
72{
73LIB"gfanlib.so";
74newstruct("divisor","ideal den,ideal num");
75newstruct("formaldivisor","list summands");
76newstruct("pdivisor","list summands, cone tail");
77
78system("install","divisor","print",divisor_print,1);
79system("install","divisor","+",divisorplus,2);
80system("install","divisor","*",proxymultdivisor,2);
81system("install","formaldivisor","print",formaldivisor_print,1);
82system("install","formaldivisor","+",formaldivisorplus,2);
83system("install","formaldivisor","*",proxymultformaldivisor,2);
84system("install","pdivisor","+",pdivisorplus,2);
85}
86
87
88
89proc divisor_print(divisor D)
90"USAGE:  divisor_print(D); D; D = divisor. @*
91ASSUME:  D is a divisor.
92RETURN:  Will print D.
93KEYWORDS: divisors
94EXAMPLE: example divisor_print; shows an example
95"
96{
97"("+string(D.num)+") - ("+string(D.den)+")";
98}
99example
100{ "EXAMPLE:";
101ring r=31991,(x,y,z),dp;
102ideal I = y^2*z - x*(x-z)*(x+3*z);
103qring Q = std(I);
104divisor P = makeDivisor(ideal(x,z),ideal(1));
105P;
106}
107
108
109proc formaldivisor_print(formaldivisor fD)
110"USAGE:  formaldivisor_print(D); D; D = formaldivisor. @*
111ASSUME:  fD is a formaldivisor.
112RETURN:  Will print fD.
113KEYWORDS: divisors
114EXAMPLE: example formaldivisor_print; shows an example
115"
116{
117  int i; string s; list L=fD.summands;
118  list cd; int c; divisor d;
119  string linebreak =
120"
121";
122  for (i=1; i<=size(L); i++)
123  {
124    cd=L[i]; c=cd[1]; d=cd[2];
125    if (i>1 && c>=0) { s = s + "+"; }
126    s = s + string(c)+"*( ("+string(d.num)+") - ("+string(d.den)+") )";
127    if (i!=size(L)) { s = s + linebreak; }
128  }
129  s;
130}
131example
132{ "EXAMPLE:";
133ring r=31991,(x,y,z),dp;
134ideal I = y^2*z - x*(x-z)*(x+3*z);
135qring Q = std(I);
136divisor P = makeDivisor(ideal(x,z),ideal(1));
137P;
138}
139
140
141////////////////////////////////////////////////////////////////////////////////
142// divisors as numerator and denomiator ideals
143
144proc makeDivisor(ideal I, ideal J)
145"USAGE:  makeDivisor(I ,J); I = ideal, J = ideal. @*
146ASSUME:  I and J are ideals in a qring Q of a smooth irreducible variety X
147         such that any ideal in Q satisfies the S2 condition.
148RETURN:  a divisor on X
149THEORY:  The procedure will eliminate all components which are not of codimension 1.
150         The S2 condition requires that every proper nonzero principal ideal
151         has pure codimension 1.
152KEYWORDS: divisors
153EXAMPLE: example makeDivisor; shows an example
154"
155{
156divisor C;
157C.num = purify1(I);
158C.den = purify1(J);
159return(C);
160}
161
162example
163{ "EXAMPLE:";
164ring r=31991,(x,y,z),dp;
165ideal I = y^2*z - x*(x-z)*(x+3*z);
166qring Q = std(I);
167divisor P = makeDivisor(ideal(x,z),ideal(1));
168}
169
170
171proc divisorplus(divisor A, divisor B)
172"USAGE:  divisorplus(A ,B); A + B; A = divisor, B = divisor. @*
173ASSUME:  A and B are divisors on X.
174RETURN:  a divisor on X
175THEORY:  The procedure will compute the product of the numerator
176         and denominator ideals, respectively.
177KEYWORDS: divisors
178EXAMPLE: example divisorplus; shows an example
179"
180{
181return(makeDivisor(interred(A.num*B.num),interred(A.den*B.den)));
182}
183example
184{ "EXAMPLE:";
185ring r=31991,(x,y,z),dp;
186ideal I = y^2*z - x*(x-z)*(x+3*z);
187qring Q = std(I);
188divisor A = makeDivisor(ideal(x,z),ideal(1));
189divisor B = makeDivisor(ideal(x,y),ideal(1));
190A+B;
191}
192
193
194proc multdivisor(int n, divisor A)
195"USAGE:  multdivisor(n ,A); A*n; n = integer, A = divisor.@*
196ASSUME:  n is an integer and A is a divisor on X.
197RETURN:  a divisor on X
198THEORY:  The procedure will compute the n-th power of the numerator
199         and denominator ideals, respectively.
200KEYWORDS: divisors
201EXAMPLE: example multdivisor; shows an example
202"
203{
204if (n<0) {return(multdivisor(-n,negativedivisor(A)));}
205return(makeDivisor(interred((A.num)^n),interred((A.den)^n)));
206}
207example
208{ "EXAMPLE:";
209ring r=31991,(x,y,z),dp;
210ideal I = y^2*z - x*(x-z)*(x+3*z);
211qring Q = std(I);
212divisor A = makeDivisor(ideal(x,z),ideal(1));
213A;
214divisor D = multdivisor(4,A);
215D;
216A*4;
217}
218
219
220/***
221 * For operator overloading, which needs a procedure which takes a divisor first
222 * and integer second.
223 **/
224proc proxymultdivisor(divisor A, int n)
225{
226if (n<0) {return(multdivisor(-n,negativedivisor(A)));}
227return(makeDivisor(interred((A.num)^n),interred((A.den)^n)));
228}
229
230
231proc negativedivisor(divisor A)
232"USAGE:  negativedivisor(A); A*(-1); A = divisor.@*
233ASSUME:  A is a divisor on X.
234RETURN:  a divisor on X
235THEORY:  The procedure will interchange the numerator and denominator ideals.
236KEYWORDS: divisors
237EXAMPLE: example negativedivisor; shows an example
238"
239{
240divisor B;
241B.num=A.den;
242B.den=A.num;
243return(B);
244}
245example
246{ "EXAMPLE:";
247ring r=31991,(x,y,z),dp;
248ideal I = y^2*z - x*(x-z)*(x+3*z);
249qring Q = std(I);
250divisor A = makeDivisor(ideal(x,z),ideal(1));
251A;
252divisor D = negativedivisor(A);
253D;
254}
255
256
257proc normalForm(divisor A)
258"USAGE:  normalForm(A); A = divisor.@*
259ASSUME:  A is a divisor on X.
260RETURN:  different representative of the same divisor on X
261THEORY:  The procedure will cancel common components of numerator and denominator.
262KEYWORDS: divisors
263EXAMPLE: example normalForm; shows an example
264"
265{
266divisor B;
267B.num=quotient(A.num,A.den);
268B.den=quotient(A.den,A.num);
269return(B);
270}
271example
272{ "EXAMPLE:";
273ring r=31991,(x,y,z),dp;
274ideal I = y^2*z - x*(x-z)*(x+3*z);
275qring Q = std(I);
276divisor A = makeDivisor(ideal(x,z),ideal(1));
277divisor B = makeDivisor(ideal(x,y),ideal(1));
278divisor D = (A+B)+multdivisor(-1,B);
279D;
280normalForm(D);
281}
282
283
284
285
286static proc isEqualIdeal(ideal A,ideal B){
287return((size(NF(A,std(B)))==0) && (size(NF(B,std(A)))==0));
288}
289
290
291proc isEqualDivisor(divisor A,divisor B)
292"USAGE:  isEqualDivisor(A,B); A = divisor, B = divisor.@*
293ASSUME:  A and B are divisors on X.
294RETURN:  int 0 or 1, checks equality of A and B.
295THEORY:  The procedure will compute the normal forms of A and B and compare.
296KEYWORDS: divisors
297EXAMPLE: example isEqualDivisor; shows an example
298"
299{
300A=normalForm(A);
301B=normalForm(B);
302return((isEqualIdeal(A.num,B.num)) && (isEqualIdeal(A.den,B.den)));
303}
304example
305{ "EXAMPLE:";
306ring r=31991,(x,y,z),dp;
307ideal I = y^2*z - x*(x-z)*(x+3*z);
308qring Q = std(I);
309divisor A = makeDivisor(ideal(x,z),ideal(1));
310divisor B = makeDivisor(ideal(x,y),ideal(1));
311divisor D = (A+B)+multdivisor(-1,B);
312isEqualDivisor(A,D);
313}
314
315
316static proc purify1(ideal I)
317{
318I = simplify(I,2);
319if (I[1]==0){ERROR("expected a non-zero ideal");}
320ideal f = I[1];
321return(minbase(quotient(f,quotient(f,I))));
322}
323
324
325static proc basis(ideal I,int d)
326{
327I=simplify(jet(intersect(I,maxideal(d)),d),2);
328return(I)}
329
330//basis(ideal(x,y^3),2);
331
332
333proc globalSections(divisor D)
334"USAGE:  globalSections(A); A = divisor.@*
335ASSUME:  A is a divisor on X.
336RETURN:  a list with a basis of the space of global sections of D.
337THEORY:  We assume that the qring of X satisfies the S2-condition and that X is smooth. We compute sat((f*J) : I) /f
338         where D = (I)-(J).
339KEYWORDS: divisors
340EXAMPLE: example globalSections; shows an example
341"
342{
343poly f =(simplify(D.num,2))[1];
344ideal LD = basis(purify1(quotient(f*D.den,D.num)),deg(f));
345list L = simplify(LD,2),f;
346return(L);
347}
348example
349{ "EXAMPLE:";
350ring r=31991,(x,y,z),dp;
351ideal I = y^2*z - x*(x-z)*(x+3*z);
352qring Q = std(I);
353divisor P = makeDivisor(ideal(x,z),ideal(1));
354divisor D = multdivisor(4,P);
355globalSections(D);
356}
357
358
359static proc sectionIdeal(poly f, poly g, divisor D){
360return(purify1(quotient(quotient(f*D.num,g), D.den)));
361}
362
363proc degreeDivisor(divisor A)
364"USAGE:  degreeDivisor(A); A = divisor.@*
365ASSUME:  A is a divisor on X.
366RETURN:  The degree of A.
367THEORY:  We compute difference of the degrees of the numerator and denominator ideals.
368KEYWORDS: divisors
369EXAMPLE: example degreeDivisor; shows an example
370"
371{
372  return( mult(std(A.num))-mult(std(A.den)));
373}
374example
375{ "EXAMPLE:";
376ring r=31991,(x,y,z),dp;
377ideal I = y^2*z - x*(x-z)*(x+3*z);
378qring Q = std(I);
379divisor P = makeDivisor(ideal(x,z),ideal(1));
380degreeDivisor(P);
381}
382
383proc linearlyEquivalent(divisor A, divisor B)
384"USAGE:  linearlyEquivalent(A,B); A = divisor, B = divisor.@*
385ASSUME:  A and B are divisors on X.
386RETURN:  list if A and B a linearly equivalent and int 0 otherwise.
387THEORY:  Checks whether A-B is principle. If yes, returns a list L=(f,g) where
388         A - B = (f/g).
389KEYWORDS: divisors
390EXAMPLE: example linearlyEquivalent; shows an example
391"
392{
393divisor F = normalForm(divisorplus(A,negativedivisor(B)));
394list LB = globalSections(F);
395if (size(LB[1])!=1) {return(0);}
396ideal V = sectionIdeal(LB[1][1,1],LB[2],F);
397if (isEqualIdeal(V,ideal(1))==1) {return(list(LB[1][1,1],LB[2]));}
398return(0);
399}
400example
401{ "EXAMPLE:";
402ring r=31991,(x,y,z),dp;
403ideal I = y^2*z - x*(x-z)*(x+3*z);
404qring Q = std(I);
405divisor A = makeDivisor(ideal(x,z),ideal(1));
406divisor B = makeDivisor(ideal(x,y),ideal(1));
407linearlyEquivalent(A,B);
408linearlyEquivalent(multdivisor(2,A),multdivisor(2,B));
409}
410
411
412proc effective(divisor A)
413"USAGE:  effective(A); A = divisor.@*
414ASSUME:  A is a divisor on X which is linearly equivalent to an effective divisor.
415RETURN:  divisor on X.
416THEORY:  We compute an effective divisor linearly equivalent to A.
417KEYWORDS: divisors
418EXAMPLE: example effective; shows an example
419"
420{
421list LB = globalSections(A);
422if (size(LB[1])==0) {ERROR("the divisor is not linearly equivalent to an effective divisor");}
423ideal V = sectionIdeal(LB[1][1,1],LB[2],A);
424return(makeDivisor(V,ideal(1)));
425}
426example
427{ "EXAMPLE:";
428ring r=31991,(x,y,z),dp;
429ideal I = y^2*z - x*(x-z)*(x+3*z);
430qring Q = std(I);
431divisor A = makeDivisor(ideal(x,z),ideal(1));
432divisor B = makeDivisor(ideal(x,y),ideal(1));
433divisor D = divisorplus(multdivisor(2,B),negativedivisor(A));
434effective(D);
435}
436
437
438
439////////////////////////////////////////////////////////////////////////////////
440// formal sums of divisors
441
442proc makeFormalDivisor(list L)
443"USAGE:  makeFormalDivisor(L); L = list.@*
444ASSUME:  L is a list of tuples of an integer and a divisor.
445RETURN:  a formal divisor on X
446THEORY:  Represents an integer formal sum of divisors.
447KEYWORDS: divisors
448EXAMPLE: example makeFormalDivisor; shows an example
449"
450{
451formaldivisor C;
452C.summands = L;
453return(C);
454}
455example
456{ "EXAMPLE:";
457ring r=31991,(x,y,z),dp;
458ideal I = y^2*z - x*(x-z)*(x+3*z);
459qring Q = std(I);
460divisor A = makeDivisor(ideal(x,z),ideal(1));
461divisor B = makeDivisor(ideal(x,y),ideal(1));
462makeFormalDivisor(list(list(-5,A),list(2,B)));
463}
464
465
466proc evaluateFormalDivisor(formaldivisor D)
467"USAGE:  evaluateFormalDivisor(D); D = formal divisor.@*
468ASSUME:  D is a formal divisor on X.
469RETURN:  a divisor on X
470THEORY:  Will evaluate the formal sum.
471KEYWORDS: divisors
472EXAMPLE: example evaluateFormalDivisor; shows an example
473"
474{
475list L = D.summands;
476if (size(L)==0) {return(makeDivisor(ideal(1),ideal(1)));}
477int i;
478divisor E = multdivisor(L[1][1],L[1][2]);
479  for ( i=2; i <= size(L); i++ )
480  {
481    E = divisorplus(E, multdivisor(L[i][1],L[i][2]));
482  }
483return(E);
484}
485example
486{ "EXAMPLE:";
487ring r=31991,(x,y,z),dp;
488ideal I = y^2*z - x*(x-z)*(x+3*z);
489qring Q = std(I);
490divisor A = makeDivisor(ideal(x,z),ideal(1));
491divisor B = makeDivisor(ideal(x,y),ideal(1));
492formaldivisor fE= makeFormalDivisor(list(list(-5,A),list(2,B)));
493evaluateFormalDivisor(fE);
494}
495
496
497
498static proc position(divisor I,list L){
499int i;
500for (i=1; i <=size(L); i++){
501   if (isEqualDivisor(I,L[i][2])==1) {return(i);}
502}
503return(0);}
504
505
506proc formaldivisorplus(formaldivisor A, formaldivisor B)
507"USAGE:  formaldivisorplus(A ,B); A + B; A = formaldivisor, B = formaldivisor. @*
508ASSUME:  A and B are formal divisors on X.
509RETURN:  a formal divisor on X
510THEORY:  The procedure will add the formal sums.
511KEYWORDS: divisors
512EXAMPLE: example formaldivisorplus; shows an example
513"
514{
515formaldivisor C;
516int i,p;
517int j=1;
518list L;
519list LA=A.summands;
520list LB=B.summands;
521for (i=1; i<=size(LA);i++){
522   p=position(LA[i][2],L);
523   if (p==0) {
524        L[j]=list();
525        L[j][2]=LA[i][2];
526        L[j][1]=LA[i][1];
527        j=j+1;
528   } else {
529        L[p][1]=L[p][1]+LA[i][1];
530   };
531}
532for (i=1; i<=size(LB);i++){
533   p=position(LB[i][2],L);
534   if (p==0) {
535        L[j]=list();
536        L[j][2]=LB[i][2];
537        L[j][1]=LB[i][1];
538        j=j+1;
539   } else {
540        L[p][1]=L[p][1]+LB[i][1];
541   };
542}
543//C.summands = (A.summands)+(B.summands);
544return(L);
545}
546example
547{ "EXAMPLE:";
548ring r=31991,(x,y,z),dp;
549ideal I = y^2*z - x*(x-z)*(x+3*z);
550qring Q = std(I);
551divisor A = makeDivisor(ideal(x,z),ideal(1));
552divisor B = makeDivisor(ideal(x,y),ideal(1));
553divisor C = makeDivisor(ideal(x-z,y),ideal(1));
554formaldivisor fE= makeFormalDivisor(list(list(-5,A),list(2,B)));
555formaldivisor fE2= makeFormalDivisor(list(list(-5,A),list(2,C)));
556formaldivisorplus(fE,fE2);
557}
558
559
560proc degreeFormalDivisor(formaldivisor A)
561"USAGE:  degreeFormalDivisor(A); A = formaldivisor.@*
562ASSUME:  A is a formaldivisor on X.
563RETURN:  The degree of A.
564THEORY:  We compute degrees of the summands and return the weighted sum.
565KEYWORDS: divisors
566EXAMPLE: example degreeFormalDivisor; shows an example
567"
568{
569int i,s;
570list L = A.summands;
571for (i=1;i<=size(L);i++){
572    s=s+L[i][1]*degreeDivisor(L[i][2]);
573}
574return(s);
575}
576example
577{ "EXAMPLE:";
578ring r=31991,(x,y,z),dp;
579ideal I = y^2*z - x*(x-z)*(x+3*z);
580qring Q = std(I);
581divisor A = makeDivisor(ideal(x,z),ideal(1));
582divisor B = makeDivisor(ideal(x,y),ideal(1));
583formaldivisor fE= makeFormalDivisor(list(list(-5,A),list(2,B)));
584degreeFormalDivisor(fE);
585}
586
587
588proc multformaldivisor(int n,formaldivisor A)
589"USAGE:  multformaldivisor(n ,A); A*n; n = integer, A = formaldivisor.@*
590ASSUME:  n is an integer and A is a formal divisor on X.
591RETURN:  a formal divisor on X
592THEORY:  The procedure will multiply the formal sum with n.
593KEYWORDS: divisors
594EXAMPLE: example multformaldivisor; shows an example
595"
596{
597formaldivisor B;
598list L=A.summands;
599int i;
600for (i=1;i<=size(L);i++){
601  L[i][1]=n*L[i][1];
602}
603B.summands=L;
604return(B);
605}
606example
607{ "EXAMPLE:";
608ring r=31991,(x,y,z),dp;
609ideal I = y^2*z - x*(x-z)*(x+3*z);
610qring Q = std(I);
611divisor A = makeDivisor(ideal(x,z),ideal(1));
612divisor B = makeDivisor(ideal(x,y),ideal(1));
613formaldivisor fE= makeFormalDivisor(list(list(-5,A),list(2,B)));
614fE*2;
615}
616
617
618/***
619 * For operator overloading, which needs a procedure which takes a divisor first
620 * and integer second.
621 **/
622proc proxymultformaldivisor(formaldivisor A, int n)
623{
624formaldivisor B;
625list L=A.summands;
626int i;
627for (i=1;i<=size(L);i++){
628  L[i][1]=n*L[i][1];
629}
630B.summands=L;
631return(B);
632}
633
634
635
636proc negativeformaldivisor(formaldivisor A)
637"USAGE:  negativeformaldivisor(A); A = formaldivisor.@*
638ASSUME:  A is a formaldivisor on X.
639RETURN:  a formal divisor on X
640THEORY:  The procedure will change the signs of the coefficients.
641KEYWORDS: divisors
642EXAMPLE: example negativeformaldivisor; shows an example
643"
644{
645formaldivisor B;
646list L=A.summands;
647int i;
648for (i=1;i<=size(L);i++){
649  L[i][1]=-L[i][1];
650}
651B.summands=L;
652return(B);
653}
654example
655{ "EXAMPLE:";
656ring r=31991,(x,y,z),dp;
657ideal I = y^2*z - x*(x-z)*(x+3*z);
658qring Q = std(I);
659divisor A = makeDivisor(ideal(x,z),ideal(1));
660divisor B = makeDivisor(ideal(x,y),ideal(1));
661formaldivisor fE= makeFormalDivisor(list(list(-5,A),list(2,B)));
662negativeformaldivisor(fE);
663}
664
665
666static proc primDecDivisor(divisor D)
667"decompose a divisor into a formal divisor of primary divisors"
668{
669formaldivisor E;
670ideal I = D.num;
671ideal J = D.den;
672list L;
673int i;
674int j = 1;
675list LI = primdecGTZ(I);
676for (i=1;i<=size(LI);i++){
677   LI[i][2];
678   L[j]=list(1,makeDivisor(LI[i][1],ideal(1)));
679   j=j+1;
680};
681list LJ = primdecGTZ(J);
682for (i=1;i<=size(LJ);i++){
683   LJ[i][2];
684   L[j]=list(-1,makeDivisor(LJ[i][1],ideal(1)));
685   j=j+1;
686};
687E.summands=L;
688return(E);}
689
690
691
692
693////////////////////////////////////////////////////////////////////////////////
694// P-divisors
695
696proc makePDivisor(list L)
697"USAGE:  makePDivisor(L); L = list.@*
698ASSUME:  L is a list of tuples of a integral polyhedron and a divisor such that
699         all polyhedra have the same tail cone.
700RETURN:  a pdivisor on X
701THEORY:  Represents an polyhedral formal sum of divisors.
702KEYWORDS: divisors, polyhedra
703EXAMPLE: example makePDivisor; shows an example
704"
705{
706pdivisor P;
707list CP = decomposePolyhedron(L[1][1]);
708P.tail = CP[1];
709list LP;
710LP[1]=list(CP[2],L[1][2]);
711int i;
712for (i=2; i<=size(L);i++){
713   CP = decomposePolyhedron(L[i][1]);
714   if (!(CP[1]==P.tail)) {ERROR("All P-coefficients should have the same tail cone");}
715   LP[i]=list(CP[2],L[i][2]);
716}
717P.summands = LP;
718return(P);
719}
720example
721{ "EXAMPLE:";
722ring r=31991,(x,y,z),dp;
723ideal I = y^2*z - x*(x-z)*(x+3*z);
724qring Q = std(I);
725divisor A = makeDivisor(ideal(x,z),ideal(1));
726divisor B = makeDivisor(ideal(x,y),ideal(1));
727intmat M[4][4]= 1,4,0,0,
728                1,0,3,0,
729                0,0,0,2,
730                1,1,1,1;
731polytope PP = polytopeViaPoints(M);
732makePDivisor(list(list(PP,A),list(PP,B)));
733}
734
735static proc decomposePolyhedron(polytope P){
736intmat rays = vertices(P);
737intmat rays2 = rays;
738int i,j;
739for (i=1; i<=nrows(rays);i++){
740   if (rays[i,1]==1) {
741      for (j=1; j<=ncols(rays);j++){
742         rays[i,j]=0;
743      }
744   } else {
745      for (j=1; j<=ncols(rays);j++){
746         rays2[i,j]=0;
747      }
748   }
749}
750cone C = coneViaPoints(rays);
751polytope C2 = polytopeViaPoints(rays2);
752return(list(C,C2));
753}
754
755
756proc evaluatePDivisor(pdivisor D,intvec v)
757"USAGE:  evaluatePDivisor(D,v); D = pdivisor, v = intvec.@*
758ASSUME:  D is a polyhedral divisor on X and v is a point in the dual of the
759         tailcone of the coefficients.
760RETURN:  a formal divisor on X
761THEORY:  Will evaluate the polyhedral sum to an integer formal sum.
762KEYWORDS: divisors, polyhedra
763EXAMPLE: example evaluatePDivisor; shows an example
764"
765{
766formaldivisor vD;
767list L = D.summands;
768cone C = D.tail;
769cone dC = dualCone(C);
770list vL;
771int i;
772intvec vv = 0,v;
773if (!(containsInSupport(dC,vv)))
774  {ERROR("the linear form given should be in the dual tail cone");}
775for (i=1; i<=size(L);i++){
776        vL[i]=list();
777        vL[i][2]=L[i][2];
778  vL[i][1]=Polymake::minimalValue(L[i][1],vv);
779}
780vD.summands = vL;
781return(vD);}
782example
783{ "EXAMPLE:";
784LIB("polymake.so");
785ring r=31991,(x,y,z),dp;
786ideal I = y^2*z - x*(x-z)*(x+3*z);
787qring Q = std(I);
788divisor A = makeDivisor(ideal(x,z),ideal(1));
789divisor B = makeDivisor(ideal(x,y),ideal(1));
790intmat M[4][4]= 1,4,0,0,
791                1,0,3,0,
792                0,0,0,2,
793                1,1,1,1;
794polytope PP = polytopeViaPoints(M);
795pdivisor pD = makePDivisor(list(list(PP,A),list(PP,B)));
796intvec v=1,1,1;
797evaluatePDivisor(pD,v);
798}
799
800
801
802
803
804proc pdivisorplus(pdivisor A, pdivisor B)
805"USAGE:  pdivisorplus(A ,B); A + B; A = pdivisor, B = pdivisor. @*
806ASSUME:  A and B are polyhedral divisors on X.
807RETURN:  a pdivisor on X
808THEORY:  The procedure will add the polyhedral formal sums by doing Minkowski sums.
809KEYWORDS: divisors, polyhedra
810EXAMPLE: example pdivisorplus; shows an example
811"
812{
813pdivisor C;
814int i,p;
815int j=1;
816if (!(A.tail==B.tail)) {ERROR("tail cones should be equal");}
817list L;
818list LA=A.summands;
819list LB=B.summands;
820for (i=1; i<=size(LA);i++){
821   p=position(LA[i][2],L);
822   if (p==0) {
823        L[j]=list();
824        L[j][2]=LA[i][2];
825        L[j][1]=LA[i][1];
826        j=j+1;
827   } else {
828        L[p][1]=Polymake::minkowskiSum(L[p][1],LA[i][1]);
829   };
830}
831for (i=1; i<=size(LB);i++){
832   p=position(LB[i][2],L);
833   if (p==0) {
834        L[j]=list();
835        L[j][2]=LB[i][2];
836        L[j][1]=LB[i][1];
837        j=j+1;
838   } else {
839        L[p][1]=Polymake::minkowskiSum(L[p][1],LB[i][1]);
840   };
841}
842C.summands = L;
843C.tail = A.tail;
844return(C);
845}
846example
847{ "EXAMPLE:";
848LIB("polymake.so");
849ring r=31991,(x,y,z),dp;
850ideal I = y^2*z - x*(x-z)*(x+3*z);
851qring Q = std(I);
852divisor A = makeDivisor(ideal(x,z),ideal(1));
853divisor B = makeDivisor(ideal(x,y),ideal(1));
854intmat M[4][4]= 1,4,0,0,
855                1,0,3,0,
856                0,0,0,2,
857                1,1,1,1;
858polytope PP = polytopeViaPoints(M);
859pdivisor pD = makePDivisor(list(list(PP,A),list(PP,B)));
860pdivisorplus(pD,pD);
861}
862
863
Note: See TracBrowser for help on using the repository browser.