source: git/Singular/LIB/scheme.lib @ c2763f

spielwiese
Last change on this file since c2763f was c2763f, checked in by Hans Schoenemann <hannes@…>, 3 years ago
fix: some fixes for exmaple isSmooth (scheme.lib)
  • Property mode set to 100644
File size: 83.2 KB
Line 
1//
2version="version scheme.lib 1.0.0.0 Apr_2021 "; // $Id$
3category="Algebraic geometry";
4info="
5LIBRARY: scheme.lib  Schemes
6
7AUTHORS:  Benjamin Mirgain
8          Janko Boehm, e-mail:boehm@mathematik.uni-kl.de
9
10OVERVIEW:
11The library implements various new classes centered around schemes.
12
13The chart type represents an affine chart of the form Spec R/I, with I being an ideal and
14R being the ring it is defined in. There is a constructor for charts which can be called by
15assigning an ideal using =.
16
17The ratFunc type represents a rational function and is made to be ring independent. This type
18is used to define morphisms between charts and schemes. You can use +, -, *, /, == for addition,
19substraction, multiplication, division and comparison. The constructor can be called by
20assigning one or two polynomials to a rational function via =.
21
22The chartmap type represents a morphism from an open subset of a chart to another chart, it
23contains two charts preim and im, representing the image and preimage of the morphism, as well as
24an ideal dom, being ideal defining the locus that is the complement of the open subset of preim
25which is the domain of the chartmap, furthermore it contains a list ratFuncs containing the
26rational functions which are the images of the variables of the ring in im under the ring morphism
27induced by the charmap. You can use * in order to compute the composition of two chartmaps, as
28long as they are compatible. The contructor of chartmaps can be called by using =.
29
30The scheme type is made up of two lists one containing the charts covering the scheme and the other
31containing the gluing morphisms between the charts.
32
33The divisor type represents a divisor of scheme.
34It contains a scheme space, representing the ambient space the divisor is in,
35as well as two lists schemes and mults, the forst containing the schemes
36which make up the divisor and the second containg the corresponding multiplicities. You can add and
37substract divisors via + and -.
38
39The morphism type represents a morphism between two schemes, it contains two schemes preim and im,
40representing the image and preimage of the morphism, as well as a list chartmaps containing the
41chartmaps between the charts of preim and the charts of im induced by the morphism. You can use *
42to compute the composition of two morphisms.
43
44Lastly the blowUpMap type is a type which represents the morphism of a blow up. In contrast two the
45morphism type it also contains all the relevant information of a blow up. It contains four schemes
46base, representing the base scheme, blow, representing the blow up, center, representing the center
47of the blow up, and exDiv, representing the exceptional divisor. Furthermore the list maps contains
48the maps contains the chartmaps defining the morphism.
49
50
51
52KEYWORDS: schemes, divisors
53
54PROCEDURES:
55blowUp(X,Y); computes the blow up of scheme X along the center Y
56ChartmapCompIdeal(phi,I); Computes the inverse image under phi of the closed subscheme defined by the ideal I
57ChartmapCompChartmap(phi,psi); Computes the composition of the chartmap phi with the chartmap psi
58ChartToScheme(U); converts chart U to a scheme
59CompleteIntersectionCover(X,W); Computes a scheme isomorphic to X, where each chart is a complete intersection
60dimension(X); computes the dimension of the scheme X
61FibreProduct(X,Y); computes the fibre product of the scheme X with the scheme Y
62HybrSmoothTest(X); smoothness test for scheme X, only works for equidimensional X
63intersection(X,Y); computes the scheme theoretic intersection of the schemes X and Y
64IntersectionNumber(D1,...); computes the intersection of the divisors D1,...
65InvImage(Y,phi); computes the inverse image of the scheme Y under the morphism phi
66IrredDec(X); computes the scheme theoretical irreducible decomposition of the scheme X
67isReduced(X); checks whether given scheme is reduced
68isSmooth(X); checks whether given scheme is smooth
69makeAffineScheme(I); constructs the affine scheme Spec R/I, where R=basering
70makeChart(I); constructor for chart type
71makeChartmap(U,V,Dom,L); constructor for chartmap type
72makeDivisor((Y1,m1),(Y2,m2),...,X); constructor for divisor type
73makeMorphism(X,Y,L); constructor for morphism type
74makeProjScheme(I); constructs the projective scheme Proj R/I, where R=basering and I homogeneous
75makeratFunc(f[,g]); constructor for ratFunc type
76MorphismCompMorphism(phi,psi); Computes the composition of the morphism phi with the morphism psi
77PolyCompratFunc(P,r_1,...,r_n); replaces the variables of the Poly P by the ratFuncs r_1,...,r_n
78ratFuncCompratFunc(F,r_1,...,r_n); replaces the variables of the ratFunc F by the ratFuncs r_1,...,r_n
79ReducedScheme(X); computes the reduced scheme of a given scheme X
80strictTransform(X,pi); strict transformation of X along the blowUpMap pi
81union(X,Y); computes the scheme theoretic union of the schemes X and Y
82";
83
84LIB "idealclass.lib";
85LIB "sing.lib";
86LIB "elim.lib";
87LIB "polyclass.lib";
88LIB "primdec.lib";
89LIB "matrix.lib";
90LIB "grobcov.lib";
91
92static proc mod_init()
93{
94  newstruct("chart","ring in, ideal value");
95  newstruct("chartmap","chart preim, chart im, ideal dom, list ratFuncs");
96  newstruct("scheme","list cover, list maps");
97  newstruct("morphism","scheme preim, scheme im, list chartmaps");
98  newstruct("divisor","scheme space, list schemes, list mults");
99  newstruct("ratFunc","ring in, poly num, poly den");
100  newstruct("blowUpMap","scheme blow, scheme base, scheme exDiv, scheme center, list maps");
101
102  system("install","chart","=",makeChart,1);
103  system("install","chartmap","=",makeChartmap,1);
104  system("install","chart","print",printChart,1);
105  system("install","chartmap","print",printChartmap,1);
106  system("install","chart","==",equalChart,2);
107  system("install","chart","!=",notequalChart,2);
108  system("install","chartmap","==",equalChartmap,2);
109  system("install","chartmap","!=",notequalChartmap,2);
110  system("install","scheme","==",equalScheme,2);
111  system("install","scheme","!=",notequalScheme,2);
112  system("install","scheme","print",printScheme,1);
113  system("install","ratFunc","=",makeratFunc,1);
114  system("install","ratFunc","print",printratFunc,1);
115  system("install","ratFunc","string",stringratFunc,4);
116  system("install","ratFunc","+",addratFunc,2);
117  system("install","ratFunc","-",subratFunc,2);
118  system("install","ratFunc","*",multratFunc,2);
119  system("install","ratFunc","/",divideratFunc,2);
120  system("install","ratFunc","==",equalratFunc,2);
121  system("install","ratFunc","!=",notequalratFunc,2);
122  system("install","ratFunc","^",expratFunc,2);
123  system("install","chartmap","[",ChartmapCompIdeal,2);
124  system("install","chartmap","*",ChartmapCompChartmap,2);
125  system("install","divisor","=",makeDivisor,1);
126  system("install","divisor","+",AddDivisor,2);
127  system("install","divisor","-",SubDivisor,2);
128  system("install","divisor","*",intTimesDivisor,2);
129  system("install","morphism","=",makeMorphism,1);
130  system("install","morphism","*",MorphismCompMorphism,2);
131}
132
133
134///////////////////////////////////////////////////
135//    Intersection Number
136///////////////////////////////////////////////////
137
138//Used for calculations in IntersectionNumber
139static proc partSum(scheme X, int i, int lower, int upper, ideal S)
140{
141    int result=0;
142    int l=lower;
143    ideal I;
144    while(l<=upper)
145    {
146          I=X.maps[i][l].dom;
147          def L=sat(S,I);
148          result=result+vdim(L[1]);
149          result=result-partSum(X,i,l+1,upper,L[1]);
150          l=l+1;
151          kill L;
152    }
153    return(result);
154}
155
156
157proc IntersectionNumber(list #)
158"USAGE: IntersectionNumber(D1,D2,...,Dn); D1,D2,... divisors
159ASSUME: The inputs all lie in general position in the same ambient space and
160        the number of inputs is equal to the dimension of the ambient space
161RETURN: The intersection number of the inputs
162EXAMPLE: example IntersectionNumber; shows an example"
163{
164    def S=basering;
165    int i,j,k,l;
166    if(typeof(#[1])!="divisor")
167    {   ERROR("Wrong inputs!");}
168    scheme X=#[1].space;
169    if(size(#)!=dimension(X))
170    {   ERROR("Wrong number of inputs!");}
171    for(k=1; k<=size(#); k=k+1)
172    {
173       if(typeof(#[k])!="divisor")
174       {   ERROR("Wrong inputs!");}
175       if(#[k].space!=X)
176       {   ERROR("Incompatible inputs!");}
177    }
178    scheme Z;
179    int result=0;
180    intvec iv;
181    int n=1;
182    for(i=1; i<=size(#); i=i+1)
183    {
184       n=n*size(#[i].schemes);
185    }
186    for(i=1;i<=n;i=i+1)
187    {
188        j=i-1;
189        for(k=1; k<=size(#); k=k+1)
190        {
191           iv[k]=(j % size(#[k].schemes))+1;
192           j=j div size(#[k].schemes);
193        }
194        j=iv[1];
195        Z=#[1].schemes[j];
196        for(k=2; k<=size(#); k=k+1)
197        {
198            j=iv[k];
199            Z=intersection(Z,#[k].schemes[j]);
200        }
201        if(dimension(Z)!=0)
202        {   ERROR("Divisors not in general position!");}
203        for(l=1;l<=size(X.cover);l=l+1)
204        {
205            def R=X.cover[l].in;
206            setring R;
207            int d=vdim(std(Z.cover[l].value+X.cover[l].value));
208            k=l+1;
209            while((d>0) and (k<=size(X.cover)))
210            {
211                ideal K=X.maps[l][k].dom;
212                def L=sat(Z.cover[l].value+X.cover[l].value,K);
213                kill K;
214                d=d-vdim(L[1]);
215                d=d+partSum(X,l,l+1,k-1,L[1]);
216                k=k+1;
217                kill L;
218            }
219            for(k=1; k<=size(#); k=k+1)
220            {
221                j=iv[k];
222                d=d*#[k].mults[j];
223            }
224            result=result+d;
225            kill d,R;
226        }
227        setring S;
228    }
229
230    setring S;
231    return(result);
232}
233example
234{
235"EXAMPLE:"; echo=2;
236ring R=0,(x,y,z),dp;
237scheme X=makeProjScheme(0);
238scheme Y=makeProjScheme(ideal(x-z),X);
239divisor A=list(Y,2),X;
240Y=makeProjScheme(ideal(x2-zy),X);
241scheme Z=makeProjScheme(ideal(x2z-z2y+y3),X);
242divisor B=list(Y,1),list(Z,1),X;
243Y=makeProjScheme(ideal(y),X);
244divisor C=list(Y,1),X;
245IntersectionNumber(3*C-A,B);
246}
247
248
249
250//////////////////////////////////////////////////////////////////////////
251//      Chart stuff
252//////////////////////////////////////////////////////////////////////////
253
254proc makeChart(ideal I)
255"USAGE: makeChart(I); ideal I
256RETURN: chart representing Spec R/I, where R=basering
257EXAMPLE: example makeChart; shows an example"
258{
259   chart C;
260   C.in=basering;
261   C.value=I;
262   return(C);
263}
264example
265{"EXAMPLE:"; echo=2;
266ring R=0,(x,y,z),dp;
267chart U=ideal(x2-yz);
268U;
269}
270
271// For comparing charts
272static proc equalChart(chart U, chart V)
273{
274   if(U.in!=V.in)
275   {
276      return(0);
277   }else{
278      def R=basering;
279      def S=U.in;
280      setring S;
281
282      if((std(quotient(U.value,V.value))!=1)&&(std(quotient(V.value,U.value))!=1))
283      {
284         setring R;
285         return(0);
286      }
287
288      setring R;
289      return(1);
290   }
291}
292
293// Also for comparing charts
294static proc notequalChart(chart U, chart V)
295{
296   return(!equalChart(U,V));
297}
298
299
300proc makeChartmap(chart Preim, chart Im, list #)
301"USAGE: makeChartmap(Preim,Im,Dom,L); Preim,Im=charts, Dom=ideal/Ideal, L=list
302        makeChartmap(Preim,Im,L); Preim,Im=charts, L=list
303ASSUME: Dom has to belong to the ring Preim.in
304        The list L is supposed to contain the entries of the list ratFuncs from the output
305        If n=size(Dom) and m=nvars(Im.in), then the list L should either be a list of n
306        lists with m ratFunc types each, or it should be a list of n*m ratFunc types in which
307        case the first m entries fill up the first list of ratFuncs and so on.
308RETURN: chartmap with the given input as data, if Dom is not given as an input, then it
309        is the ideal dom from the output will represent the largest possible domain
310EXAMPLE: example makeChartmap; creates an example"
311{
312   chartmap result;
313   result.preim=Preim;
314   result.im=Im;
315   int i,j;
316   intvec iv;
317
318   def R=basering;
319   def S=Preim.in;
320   def SS=Im.in;
321
322   setring S;
323   ideal Dom;
324   poly f;
325   if((typeof(#[1])=="ideal")||(typeof(#[1])=="Ideal"))
326   {
327      if(typeof(#[1])=="Ideal")
328      {
329         Dom=#[1].value;
330      }
331      else
332      {
333         Dom=#[1];
334      }
335      result.dom=Dom;
336      #=delete(#,1);
337      if(typeof(#[1])=="list")
338      {
339         if(typeof(#[1][1])=="ratFunc")
340         {
341            if(size(#)>1)
342            {
343            if(size(#)!=size(Dom))
344            {   ERROR("Incorrect number of inputs!");}
345            result.ratFuncs=#;
346            }else{
347            #=#[1];
348            }
349         }else{
350         #=#[1];
351         if(typeof(#[1])=="list")
352         {
353         if(size(#)!=size(Dom))
354         {   ERROR("Incorrect number of inputs!");}
355         result.ratFuncs=#;
356         }}
357      }
358      if(typeof(#[1])=="ratFunc")
359      {
360         if(size(#)!=size(Dom)*nvars(SS))
361         {   ERROR("Incorrect number of inputs!");}
362         for(i=1; i<=size(Dom); i=i+1)
363         {
364            result.ratFuncs[i]=list();
365            for(j=1; j<=nvars(SS); j=j+1)
366            {
367               result.ratFuncs[i][j]=#[(i-1)*nvars(SS)+j];
368            }
369         }
370      }
371   }
372   else
373   {
374      if(typeof(#[1])=="list")
375      {
376         f=1;
377         for(j=1; j<=size(#); j=j+1)
378         {
379            f=1;
380            for(i=1; i<=size(#[j]); i=i+1)
381            {
382                f=lcm(f,#[j][i].den);
383            }
384            Dom=Dom+f;
385         }
386         result.dom=Dom;
387         result.ratFuncs=#;
388      }
389      if(typeof(#[1])=="ratFunc")
390      {
391         if((size(#)%nvars(SS))!=0)
392         {   ERROR("Incorrect number of inputs!");}
393         int s=size(#) div nvars(SS);
394         for(i=1; i<=s; i=i+1)
395         {
396            result.ratFuncs[i]=list();
397            f=1;
398            for(j=1; j<=nvars(SS); j=j+1)
399            {
400               f=lcm(f,#[(i-1)*nvars(SS)+j].den);
401               result.ratFuncs[i][j]=#[(i-1)*nvars(SS)+j];
402            }
403            Dom=Dom+f;
404         }
405         result.dom=Dom;
406      }
407   }
408
409   kill Dom,f;
410   setring R;
411   return(result);
412}
413example
414{"EXAMPLE:"; echo=2;
415ring R=0,(a,b,c,d),dp;
416chart U=ad-bc;
417Ideal I=ideal(c,d);
418ratFunc F=a,c;
419ratFunc G=b,d;
420list L=list(list(F),list(G));
421ring S=0,(x),dp;
422chart V=0;
423
424chartmap phi=U,V,I,L;
425phi;
426
427chartmap psi=U,V,I,F,G;
428psi==phi;
429
430psi=U,V,L;
431psi==phi;
432
433psi=U,V,F,G;
434psi==phi;
435}
436
437// Print routine for chart
438static proc printChart(chart C)
439{
440   def R=basering;
441   def S=C.in;
442   setring S;
443   "ring = "+string(S);
444   "quotient = "+string(C.value);
445   setring R;
446}
447
448// Print routine for chartmap
449static proc printChartmap(chartmap phi)
450{
451   def B=basering;
452   def R=phi.preim.in;
453   def S=phi.im.in;
454   int i,k;
455
456   setring R;
457   "Preimage ring: ("+string(R)+")/("+string(phi.preim.value)+")";
458   def domR=phi.r_dom;
459   setring domR;
460   "Domain: D("+string(phi.dom)+")";
461   setring S;
462   "Image ring: ("+string(S)+")/("+string(phi.im.value)+")";
463   string st;
464   int n=size(phi.ratFuncs);
465   int m=nvars(S);
466   for(k=1; k<=n; k=k+1)
467   {
468     setring domR;
469     "Map restricted to D("+string(phi.dom[k])+"):";
470     st="     ("+varstr(R)+") -----> ( ";
471     setring R;
472     for(i=1; i<m; i=i+1)
473     {
474       st=st+string(phi.ratFuncs[k][i])+" , ";
475     }
476     st=st+string(phi.ratFuncs[k][m])+" )";
477     st;
478   }
479   setring B;
480}
481
482// Used for comparing chartmaps
483static proc equalChartmap(chartmap phi, chartmap psi)
484{
485   if((phi.preim!=psi.preim)||(phi.im!=psi.im))
486   {
487      return(0);
488   }
489
490   def R=basering;
491   def S=phi.preim.in;
492   setring S;
493
494   if(size(phi.dom)!=size(psi.dom))
495   {
496      setring R;
497      return(0);
498   }
499
500   int i,j;
501   for(i=1; i<=size(phi.dom); i=i+1)
502   {
503      if(phi.dom[i]!=psi.dom[i])
504      {
505         setring R;
506         return(0);
507      }
508   }
509
510
511   for(i=1; i<=size(phi.dom); i=i+1)
512   {
513      for(j=1; j<=size(phi.ratFuncs[i]); j=j+1)
514      {
515         poly f=phi.ratFuncs[i][j].num*psi.ratFuncs[i][j].den-psi.ratFuncs[i][j].num*phi.ratFuncs[i][j].den;
516         if(std(reduce(f,sat(phi.preim.value,ideal(phi.dom[i]))[1]))!=0)
517         {
518            setring R;
519            return(0);
520         }
521         kill f;
522      }
523   }
524
525   setring R;
526   return(1);
527}
528
529// Used for comparing chartmaps
530static proc notequalChartmap(chartmap phi, chartmap psi)
531{
532   return(!equalChartmap(phi,psi));
533}
534
535// Used for comparing schemes
536static proc equalScheme(scheme X, scheme Y)
537{
538   int i,j;
539
540   if(size(X.cover)!=size(Y.cover))
541   {
542      return(0);
543   }
544
545   for(i=1; i<=size(X.cover); i=i+1)
546   {
547      if(X.cover[i]!=Y.cover[i])
548      {
549         return(0);
550      }
551   }
552
553   for(i=1; i<=size(X.cover); i=i+1)
554   {
555      for(j=1; j<=size(X.cover); j=j+1)
556      {
557         if(X.maps[i][j]!=Y.maps[i][j])
558         {
559            return(0);
560         }
561      }
562   }
563
564   return(1);
565}
566
567// Used for comparing schemes
568proc notequalScheme(scheme X, scheme Y)
569{
570   return(!equalScheme(X,Y));
571}
572
573// Print routine for schemes
574static proc printScheme(scheme X)
575{
576   def R=basering;
577   int i,j,k,m,l;
578   for(i=1; i<=size(X.cover); i=i+1)
579   {
580      "Chart U"+string(i)+":";
581      X.cover[i];
582   }
583   "";
584   chartmap phi;
585   for(i=1; i<=size(X.cover); i=i+1)
586   {
587      for(j=1; j<=size(X.cover); j=j+1)
588      {
589         "Map U"+string(i)+" -> U"+string(j)+":";
590
591         phi=X.maps[i][j];
592         def T=phi.preim.in;
593         def S=phi.im.in;
594
595         setring T;
596         "Domain: D("+string(phi.dom)+")";
597
598         setring S;
599         string st;
600         m=nvars(S);
601         for(k=1; k<=size(phi.ratFuncs); k=k+1)
602         {
603         setring T;
604         "Map restricted to D("+string(phi.dom[k])+"):";
605         st="     ("+varstr(T)+") -----> ( ";
606         for(l=1; l<m; l=l+1)
607         {
608            st=st+string(phi.ratFuncs[k][l])+" , ";
609         }
610         st=st+string(phi.ratFuncs[k][m])+" )";
611         st;
612         "";
613         }
614         kill st;
615         kill T,S;
616      }
617   }
618   setring R;
619}
620
621
622proc ChartToScheme(chart U)
623"USAGE: ChartToScheme(U); chart U
624RETURN: scheme type corresponding to the affine scheme that corresponds to U
625EXAMPLE: example ChartToScheme; creates an example"
626{
627    def R=basering;
628    scheme result;
629
630    result.cover[1]=U;
631    def S=U.in;
632    setring S;
633
634    list L;
635    ratFunc r;
636    for(int i=1;i<=nvars(S);i=i+1)
637    {
638       r=var(i);
639       L[i]=r;
640    }
641    chartmap Phi=U,U,L;
642    result.maps=list(list(Phi));
643
644    setring R;
645    return(result);
646}
647example
648{"EXAMPLE:"; echo=2;
649ring R=0,(x,y),dp;
650chart U=x2-y3;
651scheme X=ChartToScheme(U);
652X;
653scheme Y=makeAffineScheme(x2-y3);
654X==Y;
655}
656
657
658
659
660
661///////////////////////////////////////////////////////////////////////////////////////////
662//  Basic constructions of schemes
663///////////////////////////////////////////////////////////////////////////////////////////
664
665
666proc makeAffineScheme(ideal I, list #)
667"USAGE: makeAffineScheme(I[,X]); I ideal, X scheme
668RETURN: affine Space V(I), if scheme X is given, then V(I) is created as a subscheme of X
669EXAMPLE: example makeAffineScheme; shows an example"
670{
671    scheme result;
672
673    if(size(#)==0)
674    {
675    chart U=I;
676    result.cover=list(U);
677
678    list L;
679    ratFunc r;
680    for(int i=1;i<=nvars(basering);i=i+1)
681    {
682       r=var(i);
683       L[i]=r;
684    }
685    chartmap Phi=U,U,L;
686    result.maps=list(list(Phi));
687    }
688
689    if((size(#)==1)&&(typeof(#[1])=="scheme"))
690    {
691    scheme X=#[1];
692    def R=basering;
693    def S=X.cover[1].in;
694    if((nvars(R)!=nvars(S))||(size(X.cover)!=1))
695    {   ERROR("Inputs not compatible");}
696
697    setring S;
698    ideal J=0;
699    for(int i=1;i<=nvars(S);i=i+1)
700    {
701       J=J,var(i);
702    }
703    J=delete(J,1);
704    map phi=R,J;
705
706    chart U=phi(I);
707    if(std(reduce(X.cover[1].value,std(U.value)))!=0)
708    {   ERROR("Inputs not compatible");}
709    result.cover=result.cover+list(U);
710    result.maps=X.maps;
711
712    result.maps[1][1].preim=U;
713    result.maps[1][1].im=U;
714
715    setring R;
716    }
717
718    return(result);
719}
720example
721{
722   "EXAMPLE:"; echo=2;
723   ring R=0,(x,y),dp;
724   scheme X=makeAffineScheme(ideal(0));
725   X.cover;
726   X.maps;
727   scheme Y=makeAffineScheme(ideal(x2-y3-y2),X);
728   Y.cover;
729   Y.maps;
730}
731
732
733
734
735
736proc makeProjScheme(ideal I, list #)
737"USAGE: makeProjScheme(I[,X]); I ideal, X scheme
738RETURN: projective Space V(I), if scheme X is given, then V(I) is created as a subscheme of X
739EXAMPLE: example makeProjScheme; shows an example"
740{
741    scheme result;
742    if(!homog(I))
743    {   ERROR("Ideal has to be homogeneous");}
744    def R=basering;
745
746    if(size(#)==0)
747    {
748    ring S;
749
750    int n=nvars(R);
751    int i,j;
752    chart U;
753    for(i=1;i<=n;i=i+1)
754    {
755       list rl=ringlist(R);
756       rl[2]=delete(rl[2],i);
757       rl[3]=list(list("dp",1:(n-1)),list("C",0));
758       S=ring(rl);
759
760       setring S;
761       ideal J=0;
762       for(j=1;j<=n;j=j+1)
763       {
764          if(j<i)
765          {   J=J,var(j);}
766          if(j==i)
767          {   J=J,1;}
768          if(j>i)
769          {   J=J,var(j-1);}
770       }
771       J=delete(J,1);
772       map phi=R,J;
773       kill J;
774       U=phi(I);
775       kill phi;
776
777       result.cover[i]=U;
778    }
779    setring R;
780
781    chartmap phi;
782    chart Preim;
783    chart Im;
784    int k;
785    list L;
786    for(i=1;i<=n;i=i+1)
787    {
788       result.maps[i]=list();
789       for(j=1;j<=n;j=j+1)
790       {
791          Preim=result.cover[i];
792          Im=result.cover[j];
793
794          S=Preim.in;
795          setring S;
796          L=list();
797          for(k=1;k<=n;k=k+1)
798          {
799             if(k!=j)
800             {
801             if(i<j)
802             {
803                if(k<i)
804                {   ratFunc rat=var(k),var(j-1);}
805                if(k==i)
806                {   ratFunc rat=1,var(j-1);}
807                if((k>i))
808                {   ratFunc rat=var(k-1),var(j-1);}
809             }
810             if(i==j)
811             {
812                if(k<j)
813                {   ratFunc rat=var(k);}
814                if(k>j)
815                {   ratFunc rat=var(k-1);}
816             }
817             if(i>j)
818             {
819                if((k<i))
820                {   ratFunc rat=var(k),var(j);}
821                if(k==i)
822                {   ratFunc rat=1,var(j);}
823                if((k>i))
824                {   ratFunc rat=var(k-1),var(j);}
825             }
826             L=L+list(rat);
827             kill rat;
828             }
829          }
830          phi=Preim,Im,L;
831
832          result.maps[i][j]=phi;
833       }
834    }
835    }
836
837    if((size(#)==1)&&(typeof(#[1])=="scheme"))
838    {
839    scheme X=#[1];
840    if(nvars(R)!=size(X.cover))
841    {   ERROR("Inputs are incompatible");}
842
843    ring S;
844    int n=nvars(R);
845    int i,j;
846    for(i=1;i<=n;i=i+1)
847    {
848       S=X.cover[i].in;
849       if(nvars(S)!=n-1)
850       {   ERROR("Inputs are incompatible");}
851
852       setring S;
853       ideal J=0;
854       for(j=1;j<=n;j=j+1)
855       {
856          if(j<i)
857          {   J=J,var(j);}
858          if(j==i)
859          {   J=J,1;}
860          if(j>i)
861          {   J=J,var(j-1);}
862       }
863       J=delete(J,1);
864       map phi=R,J;
865       kill J;
866
867       chart U=phi(I);
868       if(reduce(X.cover[i].value,std(U.value))!=0)
869       {   ERROR("Inputs are incompatible");}
870       result.cover[i]=U;
871       kill phi,U;
872    }
873
874    result.maps=X.maps;
875    for(i=1;i<=n;i=i+1)
876    {
877       for(j=1;j<=n;j=j+1)
878       {
879          result.maps[i][j].preim=result.cover[i];
880          result.maps[i][j].im=result.cover[j];
881       }
882    }
883    }
884
885    setring R;
886    return(result);
887}
888example
889{
890   "EXAMPLE:"; echo=2;
891   ring R=0,(x,y,z),dp;
892   scheme X=makeProjScheme(ideal(0));
893   X;
894   scheme Y=makeProjScheme(ideal(x2z-y3-y2z),X);
895   Y;
896}
897
898
899
900proc FibreProduct(scheme X, scheme Y)
901"USAGE: FibreProduct(X,Y); X,Y scheme
902RETURN: The Fibre Product of X and Y
903EXAMPLE: example FibreProduct; shows an example"
904{
905   scheme result;
906   def R=basering;
907   int i,j,n,m,k,l;
908   intvec iv,jv;
909
910   for(i=1; i<=size(X.cover); i=i+1)
911   {
912   for(j=1; j<=size(Y.cover); j=j+1)
913   {
914      def Sx=X.cover[i].in;
915      def Sy=Y.cover[j].in;
916      def S=ringtensor(Sx,Sy);
917      setring S;
918
919      n=nvars(Sx);
920      m=nvars(Sy);
921      ideal Xvars,Yvars;
922      for(k=1; k<=n; k=k+1)
923      {
924         Xvars=Xvars+var(k);
925      }
926      for(k=1; k<=m; k=k+1)
927      {
928         Yvars=Yvars+var(n+k);
929      }
930      map fx=Sx,Xvars;
931      map fy=Sy,Yvars;
932
933      setring Sx;
934      ideal Ix=X.cover[i].value;
935      setring Sy;
936      ideal Iy=Y.cover[j].value;
937      setring S;
938      chart U=fx(Ix)+fy(Iy);
939      result.cover=result.cover+list(U);
940
941      setring Sx;
942      kill Ix;
943      setring Sy;
944      kill Iy;
945      setring S;
946      kill U,fx,fy,Xvars,Yvars;
947      kill Sx,Sy,S;
948   }}
949
950   list Dx,Dy,Xpartmaps,Ypartmaps,Partmaps,Maps;
951   chartmap Phi;
952   for(i=1; i<=size(X.cover); i=i+1)
953   {
954   for(j=1; j<=size(Y.cover); j=j+1)
955   {
956   Maps=list();
957   for(k=1; k<=size(X.cover); k=k+1)
958   {
959   for(l=1; l<=size(Y.cover); l=l+1)
960   {
961      Dx=list();
962      Dy=list();
963      Xpartmaps=list();
964      Ypartmaps=list();
965      Partmaps=list();
966
967      def Sx=X.cover[i].in;
968      def Sy=Y.cover[j].in;
969      def S=result.cover[(i-1)*size(Y.cover)+j].in;
970      setring S;
971
972      ideal Xvars,Yvars,Dom;
973      for(n=1; n<=nvars(Sx); n=n+1)
974      {
975         Xvars=Xvars+var(n);
976      }
977      for(n=1; n<=nvars(Sy); n=n+1)
978      {
979         Yvars=Yvars+var(n+nvars(Sx));
980      }
981      map fx=Sx,Xvars;
982      map fy=Sy,Yvars;
983
984      Poly P;
985      ratFunc F;
986      for(n=1; n<=size(X.maps[i][k].ratFuncs); n=n+1)
987      {
988         Xpartmaps[n]=list();
989         setring Sx;
990         poly f,g;
991         setring S;
992         for(m=1; m<=size(X.maps[i][k].ratFuncs[n]); m=m+1)
993         {
994            setring Sx;
995            f=X.maps[i][k].ratFuncs[n][m].num;
996            g=X.maps[i][k].ratFuncs[n][m].den;
997            setring S;
998            F=fx(f),fx(g);
999            Xpartmaps[n][m]=F;
1000         }
1001         setring Sx;
1002         f=X.maps[i][k].dom[n];
1003         setring S;
1004         Dx[n]=fx(f);
1005
1006         setring Sx;
1007         kill f,g;
1008         setring S;
1009      }
1010
1011      for(n=1; n<=size(Y.maps[j][l].ratFuncs); n=n+1)
1012      {
1013         Ypartmaps[n]=list();
1014         setring Sy;
1015         poly f,g;
1016         setring S;
1017         for(m=1; m<=size(Y.maps[j][l].ratFuncs[n]); m=m+1)
1018         {
1019            setring Sy;
1020            f=Y.maps[j][l].ratFuncs[n][m].num;
1021            g=Y.maps[j][l].ratFuncs[n][m].den;
1022            setring S;
1023            F=fy(f),fy(g);
1024            Ypartmaps[n][m]=F;
1025         }
1026         setring Sy;
1027         f=Y.maps[j][l].dom[n];
1028         setring S;
1029         Dy[n]=fy(f);
1030
1031         setring Sy;
1032         kill f,g;
1033         setring S;
1034      }
1035
1036      for(n=1; n<=size(X.maps[i][k].ratFuncs); n=n+1)
1037      {
1038      for(m=1; m<=size(Y.maps[j][l].ratFuncs); m=m+1)
1039      {
1040         Dom=Dom+Dx[n]*Dy[m];
1041         Partmaps=Partmaps+Xpartmaps[n]+Ypartmaps[m];
1042      }}
1043
1044      Phi=result.cover[(i-1)*size(Y.cover)+j],result.cover[(k-1)*size(Y.cover)+l],Dom,Partmaps;
1045      Maps=Maps+list(Phi);
1046
1047      setring S;
1048      kill P,F,fx,fy,Xvars,Yvars,Dom;
1049      kill Sx,Sy,S;
1050   }}
1051   result.maps=result.maps+list(Maps);
1052   }}
1053
1054   return(result);
1055}
1056example
1057{"EXAMPLE:"; echo=2;
1058ring R=0,(x(1..3)),dp;
1059scheme X=makeProjScheme(ideal(x(1)^2-x(2)*x(3)));
1060ring S=0,(y(1..2)),dp;
1061scheme Y=makeProjScheme(ideal(y(1)-y(2)));
1062
1063scheme Z=FibreProduct(X,Y);
1064Z;
1065}
1066
1067
1068proc intersection(scheme X, scheme Y)
1069"USAGE: intersection(X,Y); X,Y scheme
1070RETURN: Scheme theoretic intersection of X and Y
1071EXAMPLE: example intersection; shows an example"
1072{
1073   if(size(X.cover)!=size(Y.cover))
1074   {   ERROR("Inputs not compatible");}
1075
1076   def R=basering;
1077   scheme result;
1078   int i,j;
1079   for(i=1; i<=size(X.cover); i=i+1)
1080   {
1081      if(X.cover[i].in!=Y.cover[i].in)
1082      {   ERROR("Inputs not compatible");}
1083
1084      def S=X.cover[i].in;
1085      setring S;
1086      chart U=X.cover[i].value+Y.cover[i].value;
1087      result.cover[i]=U;
1088
1089      kill U;
1090      kill S;
1091   }
1092
1093   result.maps=X.maps;
1094
1095   for(i=1; i<=size(result.cover); i=i+1)
1096   {
1097      for(j=1; j<=size(result.cover); j=j+1)
1098      {
1099         result.maps[i][j].preim=result.cover[i];
1100         result.maps[i][j].im=result.cover[j];
1101      }
1102   }
1103
1104   setring R;
1105   return(result);
1106}
1107example
1108{"EXAMPLE:"; echo=2;
1109ring R=0,(x,y,z),dp;
1110scheme X=makeProjScheme(0);
1111scheme Y=makeProjScheme(ideal(x2-yz),X);
1112scheme Z=makeProjScheme(ideal(xy2-z3),X);
1113
1114scheme S=intersection(Y,Z);
1115S;
1116}
1117
1118
1119
1120proc union(scheme X, scheme Y)
1121"USAGE: union(X,Y); X,Y scheme
1122RETURN: Scheme theoretic union of X and Y
1123EXAMPLE: example union; shows an example"
1124{
1125   if(size(X.cover)!=size(Y.cover))
1126   {   ERROR("Inputs not compatible");}
1127
1128   def R=basering;
1129   scheme result;
1130   int i,j;
1131   for(i=1; i<=size(X.cover); i=i+1)
1132   {
1133      if(X.cover[i].in!=Y.cover[i].in)
1134      {   ERROR("Inputs not compatible");}
1135
1136      def S=X.cover[i].in;
1137      setring S;
1138      chart U=intersect(X.cover[i].value,Y.cover[i].value);
1139      result.cover[i]=U;
1140
1141      kill U;
1142      kill S;
1143   }
1144
1145   result.maps=X.maps;
1146
1147   for(i=1; i<=size(result.cover); i=i+1)
1148   {
1149      for(j=1; j<=size(result.cover); j=j+1)
1150      {
1151         result.maps[i][j].preim=result.cover[i];
1152         result.maps[i][j].im=result.cover[j];
1153      }
1154   }
1155
1156   setring R;
1157   return(result);
1158}
1159example
1160{"EXAMPLE:"; echo=2;
1161ring R=0,(x,y,z),dp;
1162scheme X=makeProjScheme(0);
1163scheme Y=makeProjScheme(ideal(x2-yz),X);
1164scheme Z=makeProjScheme(ideal(xy2-z3),X);
1165
1166scheme S=union(Y,Z);
1167S;
1168}
1169
1170
1171proc ReducedScheme(scheme X)
1172"USAGE: ReducedScheme(X); X scheme
1173RETURN: The associated reduced scheme
1174EXAMPLE: example ReducedScheme; shows an example"
1175{
1176   def R=basering;
1177   scheme result;
1178   int i,j;
1179   for(i=1; i<=size(X.cover); i=i+1)
1180   {
1181      def S=X.cover[i].in;
1182      setring S;
1183      chart U=radical(X.cover[i].value);
1184      result.cover[i]=U;
1185
1186      kill U;
1187      kill S;
1188   }
1189
1190   result.maps=X.maps;
1191
1192   for(i=1; i<=size(result.cover); i=i+1)
1193   {
1194      for(j=1; j<=size(result.cover); j=j+1)
1195      {
1196         result.maps[i][j].preim=result.cover[i];
1197         result.maps[i][j].im=result.cover[j];
1198      }
1199   }
1200
1201   setring R;
1202   return(result);
1203}
1204example
1205{"EXAMPLE:"; echo=2;
1206ring R=0,(x,y,z),dp;
1207scheme X=makeProjScheme(x2y3z4);
1208"Original Scheme:";
1209X;
1210
1211"Reduced Scheme:";
1212scheme S=ReducedScheme(X);
1213S;
1214}
1215
1216
1217
1218
1219
1220
1221
1222
1223//////////////////////////////////////////////////////////////////////////////////////////
1224//  Schemes Miscelanious
1225//////////////////////////////////////////////////////////////////////////////////////////
1226
1227proc dimension(scheme X)
1228"USAGE: dimension(X); X scheme
1229RETURN: The dimension of the scheme X
1230EXAMPLE: example dimension; shows an example"
1231{
1232    def R=basering;
1233    int d=-1;
1234    for(int i=1; i<=size(X.cover); i=i+1)
1235    {
1236       def S=X.cover[i].in;
1237       setring S;
1238       if(d<dim(std(X.cover[i].value)))
1239       {   d=dim(std(X.cover[i].value));}
1240       kill S;
1241    }
1242
1243    setring R;
1244    return(d);
1245}
1246example
1247{"EXAMPLE:"; echo=2;
1248ring R=0,(x,y,z),dp;
1249scheme X=makeAffineScheme(x-yz);
1250dimension(X);
1251
1252X=makeProjScheme(x);
1253dimension(X);
1254
1255X=makeProjScheme(ideal(x,y));
1256dimension(X);
1257
1258X=makeProjScheme(ideal(x,y,z));
1259dimension(X);
1260}
1261
1262
1263
1264proc IrredDec(scheme X)
1265"USAGE: IrredDec(X); X scheme
1266RETURN: A list L of the irreducible components and their associated points
1267            L[i][1] is the i-th irreducible component
1268            L[i][2] is the associated point of the i-th component
1269EXAMPLE: example IrredDec; shows an example"
1270{
1271    def R=basering;
1272    list Dec;
1273
1274    int i,j,k,t;
1275    for(i=1; i<=size(X.cover); i=i+1)
1276    {
1277       def S=X.cover[i].in;
1278       setring S;
1279       list L=primdecGTZ(X.cover[i].value);
1280       for(j=1; j<=size(L); j=j+1)
1281       {
1282          Ideal Q,P;
1283          Q=L[j][1];
1284          P=L[j][2];
1285          scheme Y,Z;
1286          if(std(Q.value)!=1)
1287          {
1288          Y.maps=X.maps;
1289          Z.maps=X.maps;
1290          t=1;
1291          for(k=1; k<=size(Dec); k=k+1)
1292          {
1293             if((std(quotient(P.value,Dec[k][2].cover[i].value))==1)&&(std(quotient(Dec[k][2].cover[i].value,P.value))==1))
1294             {
1295                t=0;
1296             }
1297          }
1298          if(t)
1299          {
1300             for(k=1; k<=size(X.cover); k=k+1)
1301             {
1302                Ideal QQ,PP;
1303                chart CQ,CP;
1304                def T=X.cover[k].in;
1305                setring T;
1306                QQ=ChartmapCompIdeal(X.maps[k][i],Q);
1307                PP=ChartmapCompIdeal(X.maps[k][i],P);
1308                CQ=QQ.value;
1309                CP=PP.value;
1310                Y.cover[k]=CQ;
1311                Z.cover[k]=CP;
1312                kill T;
1313                kill CQ,CP,QQ,PP;
1314             }
1315             Dec=Dec+list(list(Y,Z));
1316          }
1317          setring S;
1318          }
1319          kill Q,P;
1320          kill Y,Z;
1321       }
1322
1323       kill L;
1324       kill S;
1325    }
1326
1327    if(size(Dec)==0)
1328    {
1329       Dec=list(list(X,X));
1330    }
1331
1332    setring R;
1333    return(Dec);
1334}
1335example
1336{"EXAMPLE:"; echo=2;
1337ring R=0,(x,y,z),dp;
1338scheme X=makeProjScheme(ideal(x2,xy));
1339list L=IrredDec(X);
1340
1341L[1][1].cover;
1342
1343L[2][1].cover;
1344
1345X==union(L[1][1],L[2][1]);
1346}
1347
1348
1349proc isReduced(scheme X)
1350"USAGE: isReduced(X); X scheme
1351RETURN: 1 if the scheme is reduced; 0 otherwise
1352EXAMPLE: example isReduced; shows an example"
1353{
1354   def R=basering;
1355   for(int i=1; i<=size(X.cover); i=i+1)
1356   {
1357      def S=X.cover[i].in;
1358      setring S;
1359      if(std(quotient(X.cover[i].value,radical(X.cover[i].value)))!=1)
1360      {
1361         setring R;
1362         return(0);
1363      }
1364
1365      kill S;
1366   }
1367
1368   setring R;
1369   return(1);
1370}
1371example
1372{"EXAMPLE:"; echo=2;
1373ring R=0,(x,y,z),dp;
1374scheme X=makeProjScheme(x2y3z4);
1375"Original Scheme is not reduced:";
1376isReduced(X);
1377
1378"Reduced Scheme is reduced:";
1379scheme S=ReducedScheme(X);
1380isReduced(S);
1381}
1382
1383
1384
1385
1386
1387
1388
1389
1390////////////////////////////////////////////////////////////////////////////////
1391//    divisor stuff
1392////////////////////////////////////////////////////////////////////////////////
1393
1394proc makeDivisor(list #)
1395"USAGE: makeDivisor((Y1,m1)[,(Y2,m2),...],X); Y1,Y2,.. scheme , m1,m2,.. int , X scheme
1396RETURN: The divisor D=m1*[Y1]+m2*[Y2]+...
1397EXAMPLE: example makeDivisor; shows an example"
1398{
1399   def R=basering;
1400   if(typeof(#[size(#)])!="scheme")
1401   { ERROR("Last entry must a scheme")}
1402   scheme X=#[size(#)];
1403
1404   divisor result;
1405   result.space=X;
1406   for(int i=1; i<=size(#)-1; i=i+1)
1407   {
1408      result.schemes[i]=#[i][1];
1409      result.mults[i]=#[i][2];
1410   }
1411
1412   setring R;
1413   return(result);
1414}
1415example
1416{"EXAMPLE:"; echo=2;
1417ring R=0,(x,y,z),dp;
1418scheme X=makeProjScheme(0);
1419
1420scheme Y=makeProjScheme(ideal(x-z),X);
1421divisor A=list(Y,2),X;
1422
1423Y=makeProjScheme(ideal(x2-zy),X);
1424scheme Z=makeProjScheme(ideal(x2z-z2y+y3),X);
1425divisor B=list(Y,1),list(Z,1),X;
1426
1427Y=makeProjScheme(ideal(y),X);
1428divisor C=list(Y,1),X;
1429
1430divisor D=A+3*B-C;
1431D.schemes;
1432D.mults;
1433
1434divisor E=2*A-4*B;
1435E.schemes;
1436E.mults;
1437
1438divisor F=D+E;
1439F.schemes;
1440F.mults;
1441}
1442
1443
1444// Addition of divisors
1445static proc AddDivisor(divisor C, divisor D)
1446{
1447    if(C.space!=D.space)
1448    {   ERROR("Divisors incompatible.");}
1449
1450    def R=basering;
1451    divisor result=C;
1452    result.space=C.space;
1453    int i,k,check;
1454    for(i=1; i<=size(D.schemes); i=i+1)
1455    {
1456       check=0;
1457       k=1;
1458       while((k<=size(result.schemes))&&(check==0))
1459       {
1460          if(D.schemes[i]==result.schemes[k])
1461          {
1462             check=1;
1463             result.mults[k]=result.mults[k]+D.mults[i];
1464             if(result.mults[k]==0)
1465             {
1466                result.schemes=delete(result.schemes,k);
1467                result.mults=delete(result.mults,k);
1468             }
1469          }
1470          k=k+1;
1471       }
1472
1473       if(check==0)
1474       {
1475          result.schemes=result.schemes+list(D.schemes[i]);
1476          result.mults=result.mults+list(D.mults[i]);
1477       }
1478    }
1479
1480    setring R;
1481    return(result);
1482}
1483
1484
1485// Substraction of divisors
1486static proc SubDivisor(divisor C, divisor D)
1487{
1488    divisor result=(-1)*D;
1489    result=result+C;
1490    return(result);
1491}
1492
1493
1494// Multiplication of a divisor with an integer
1495static proc intTimesDivisor(int n, divisor D)
1496{
1497    divisor result;
1498    result.space=D.space;
1499    if(n!=0)
1500    {
1501    result.schemes=D.schemes;
1502    for(int i=1; i<=size(D.schemes); i=i+1)
1503    {
1504       result.mults[i]=n*D.mults[i];
1505    }}
1506    return(result);
1507}
1508
1509
1510
1511
1512
1513
1514
1515////////////////////////////////////////////////////////////////////////////////
1516//    ratFunc stuff
1517////////////////////////////////////////////////////////////////////////////////
1518
1519proc makeratFunc(poly num, list #)
1520"USAGE: makeratFunc(f[,g]); poly f,g
1521RETURN: rational function f/g, if no g is put in then g is set to 1 by default
1522EXAMPLE: example makeratFunc; creates an example"
1523{
1524ratFunc r;
1525r.in=basering;
1526if(size(#)==0)
1527{
1528   r.num=num;
1529   r.den=1;
1530} else {
1531   if(#[1]==0)
1532   {   ERROR("Denominator cannot be 0!");}
1533   poly g=gcd(num,#[1]);
1534   g=g*leadcoef(#[1]);
1535   r.num=num/g;
1536   r.den=#[1]/g;
1537}
1538return(r);
1539}
1540example
1541{"EXAMPLE:"; echo=2;
1542ring R=0,(x,y),dp;
1543ratFunc F=x2,xy;
1544F;
1545ratFunc G=x2+y3;
1546G;
1547ratFunc H=1,x;
1548H;
1549
1550G+H-F;
1551F*G;
1552F^3;
1553F/G;
1554
1555F==H;
15561/(y*F)==H;
1557}
1558
1559// Print routine for ratFunc
1560static proc printratFunc(ratFunc r)
1561{
1562   def R=basering;
1563   def SS=r.in;
1564   setring SS;
1565   if(r.den==1)
1566   {
1567      string(r.num);
1568   } else {
1569      string(r.num)+" / "+string(r.den);
1570   }
1571   setring R;
1572}
1573
1574
1575// Gives ratFunc as a string
1576static proc stringratFunc(ratFunc r)
1577{
1578   def R=basering;
1579   def SS=r.in;
1580   setring SS;
1581   if(r.den==1)
1582   {
1583      string st=string(r.num);
1584   } else {
1585      string st=string(r.num)+" / "+string(r.den);
1586   }
1587   setring R;
1588   return(st);
1589}
1590
1591
1592// Used for comparing ratFuncs
1593static proc equalratFunc(ratFunc r, ratFunc q)
1594{
1595   def R=basering;
1596   def SS=r.in;
1597   setring SS;
1598   int requalq = (r.num*q.den)==(r.den*q.num);
1599   setring R;
1600   return(requalq);
1601}
1602
1603
1604// Used for comparing ratFuncs
1605static proc notequalratFunc(ratFunc F, ratFunc G)
1606{
1607   return(!equalratFunc(F,G));
1608}
1609
1610
1611// Multiplication of ratFuncs
1612static proc multratFunc(ratFunc r, ratFunc q)
1613{
1614   def R=basering;
1615   def SS=r.in;
1616   setring SS;
1617   ratFunc rtimesq = r.num*q.num, r.den*q.den;
1618   setring R;
1619   return(rtimesq);
1620}
1621
1622
1623// Division of ratFuncs
1624static proc divideratFunc(ratFunc r, ratFunc q)
1625{
1626   def R=basering;
1627   def SS=r.in;
1628   setring SS;
1629   ratFunc zero=0,1;
1630   if(q==zero)
1631   {  ERROR("Cannot divide by zero!!!")}
1632   ratFunc rdivq = r.num*q.den, r.den*q.num;
1633   setring R;
1634   return(rdivq);
1635}
1636
1637
1638// Addition of ratFuncs
1639static proc addratFunc(ratFunc r, ratFunc q)
1640{
1641   def R=basering;
1642   def SS=r.in;
1643   setring SS;
1644   ratFunc rplusq = r.num*q.den+r.den*q.num, r.den*q.den;
1645   setring R;
1646   return(rplusq);
1647}
1648
1649
1650// Substraction of ratFuncs
1651static proc subratFunc(ratFunc r, ratFunc q)
1652{
1653   def R=basering;
1654   def SS=r.in;
1655   setring SS;
1656   ratFunc rminusq = r.num*q.den-r.den*q.num, r.den*q.den;
1657   setring R;
1658   return(rminusq);
1659}
1660
1661
1662// Exponents of ratFuncs
1663static proc expratFunc(ratFunc r, int n)
1664{
1665   def R=basering;
1666   def SS=r.in;
1667   setring SS;
1668   ratFunc zero=0,1;
1669   ratFunc rexpn;
1670   if(n>=0)
1671   {
1672      rexpn=r.num^n,r.den^n;
1673   } else {
1674      if(r==zero)
1675      {   ERROR("Negative exponent not allowed for 0!!!");}
1676      rexpn=r.den^(-n),r.num^(-n);
1677   }
1678   setring R;
1679   return(rexpn);
1680}
1681
1682
1683
1684
1685////////////////////////////////////////////////////////////////////////////////
1686//   Composition stuff
1687////////////////////////////////////////////////////////////////////////////////
1688
1689
1690proc PolyCompratFunc(Poly F, list #)
1691"USAGE: PolyCompratFunc(F,r_1,...,r_n); F Poly, r_1,...,r_n ratFunc
1692ASSUME: r_1,...,r_n belong to the same ring and n=nvars(F.in)
1693RETURN: F(r_1,...,r_n)
1694EXAMPLE: example PolyCompratFunc; creates an example"
1695{
1696   def R=basering;
1697   def S=F.in;
1698   setring S;
1699   poly Fp=F.value;
1700   int n=size(#);
1701   if(n!=nvars(S))
1702   {   ERROR("Wrong number of arguments!!!");}
1703
1704   list rl=ringlist(S);
1705   list newvars=list();
1706   for(int k=1; k<=n; k=k+1)
1707   {
1708      newvars=newvars+list("x("+string(k)+")");
1709   }
1710   newvars=newvars+list("t");
1711   rl[2]=newvars;
1712   rl[3][1][2]=rl[3][1][2],1;
1713   def newS=ring(rl);
1714   setring newS;
1715
1716   ideal I=var(1);
1717   for(k=2; k<=n; k=k+1)
1718   {
1719      I=I,ideal(var(k));
1720   }
1721   map phi=S,I;
1722   poly f=phi(Fp);
1723   f=homog(f,t);
1724   int d=deg(f);
1725
1726   def SS=#[1].in;
1727   setring SS;
1728   poly g=1;
1729   for(k=1; k<=n; k=k+1)
1730   {
1731      g=g*#[k].den;
1732   }
1733   ideal J=(g/#[1].den)*#[1].num;
1734   for(k=2; k<=n; k=k+1)
1735   {
1736      J=J,ideal((g/#[k].den)*#[k].num);
1737   }
1738   J=J,ideal(g);
1739   map psi=newS,J;
1740   ratFunc result=psi(f),g^d;
1741
1742   setring R;
1743   return(result);
1744}
1745example
1746{"EXAMPLE:"; echo=2;
1747ring R=0,(x,y),dp;
1748Poly F=x2-xy+y3;
1749
1750ring S=0,(u,v,w),dp;
1751ratFunc r1=u,v;
1752ratFunc r2=w,u;
1753
1754ratFunc G=PolyCompratFunc(F,r1,r2);
1755G;
1756}
1757
1758
1759proc ratFuncCompratFunc(ratFunc r, list #)
1760"USAGE: ratFuncCompratFunc(F,r_1,...,r_n); F ratFunc, r_1,...,r_n ratFunc
1761ASSUME: r_1,...,r_n belong to the same ring and n=nvars(F.in)
1762RETURN: F(r_1,...,r_n)
1763EXAMPLE: example ratFuncCompratFunc; creates an example"
1764{
1765   def R=basering;
1766   def S=r.in;
1767   setring S;
1768   if(size(#)!=nvars(S))
1769   {   ERROR("Wrong number of arguments!!!");}
1770   ratFunc result=PolyCompratFunc(r.num,#)/PolyCompratFunc(r.den,#);
1771   setring R;
1772   return(result);
1773}
1774example
1775{"EXAMPLE:"; echo=2;
1776ring R=0,(x,y),dp;
1777ratFunc F=x2+y,x-y;
1778
1779ring S=0,(u,v,w),dp;
1780ratFunc r1=u,v;
1781ratFunc r2=u+w,v-w;
1782
1783ratFunc G=ratFuncCompratFunc(F,r1,r2);
1784G;
1785}
1786
1787
1788
1789proc ChartmapCompIdeal(chartmap phi, Ideal I)
1790"USAGE: ChartmapCompIdeal(phi,I); phi chartmap, I ideal
1791ASSUME: I.in=phi.im.in
1792RETURN: Ideal defining the inverse image under phi of the closed subscheme defined by I
1793EXAMPLE: example ChartmapCompIdeal; creates an example"
1794{
1795   def R=basering;
1796   def S=phi.preim.in;
1797   def SS=phi.im.in;
1798   setring SS;
1799   int m=size(phi.ratFuncs);
1800   int n=size(I.value);
1801   int i,j;
1802   ratFunc r;
1803   list L,LL;
1804   Ideal II;
1805   def domR=phi.r_dom;
1806   for(i=1; i<=m; i=i+1)
1807   {
1808      setring SS;
1809      LL=list();
1810      for(j=1;j<=n;j=j+1)
1811      {
1812         r=PolyCompratFunc(I.value[j],phi.ratFuncs[i]);
1813         LL=LL+list(r);
1814      }
1815      setring S;
1816      poly k=0;
1817      for(j=1;j<=n;j=j+1)
1818      {
1819         k=lcm(k,LL[j].den);
1820      }
1821      II=0;
1822      poly f;
1823      for(j=1;j<=n;j=j+1)
1824      {
1825         f=k/LL[j].den;
1826         II=II+ideal(f*LL[j].num);
1827      }
1828      kill f,k;
1829      II.value=II.value+phi.preim.value;
1830      setring domR;
1831      poly tmp=phi.dom[i];
1832      setring S;
1833      if (S!=domR) { poly tmp=fetch(domR,tmp);}
1834      II=sat(II.value,tmp)[1];
1835      kill tmp;
1836      L[i]=II;
1837      if (S!=domR) { setring domR;kill tmp;setring S;}
1838   }
1839
1840   setring S;
1841   ideal J=1;
1842   Ideal result;
1843   for(i=1; i<=m; i=i+1)
1844   {
1845      J=intersect(J,L[i].value);
1846   }
1847   result=std(J);
1848
1849   setring R;
1850   return(result);
1851}
1852example
1853{
1854ring R=0,(a,b,c,d),dp;
1855chart U=ideal(ad-bc,ac-bd);
1856Ideal Dom=ideal(c,d);
1857ratFunc F1=a,c;
1858ratFunc F2=b,c;
1859ratFunc G1=b,d;
1860ratFunc G2=a,d;
1861ring S=0,(x,y),dp;
1862chart V=0;
1863Ideal I=ideal(x,y);
1864chartmap phi=U,V,Dom,F1,F2,G1,G2;
1865
1866Ideal J=ChartmapCompIdeal(phi,I);
1867J;
1868}
1869
1870
1871
1872proc ChartmapCompChartmap(chartmap phi, chartmap psi)
1873{
1874   if(phi.preim.in!=psi.im.in)
1875   {   ERROR("Chartmaps not compatible");}
1876
1877   def R=basering;
1878   def S=psi.preim.in;
1879   def SS=phi.preim.in;
1880   list Maps;
1881   int i,j,k;
1882   ratFunc G,F;
1883   setring S;
1884   ideal Dom;
1885   int l=1;
1886   for(i=1; i<=size(psi.ratFuncs); i=i+1)
1887   {
1888      for(j=1; j<=size(phi.ratFuncs); j=j+1)
1889      {
1890         setring SS;
1891         G=phi.dom[j];
1892         F=ratFuncCompratFunc(G,psi.ratFuncs[i]);
1893         setring S;
1894         Dom[l]=lcm(F.num,psi.dom[i]);
1895         for(k=1; k<=nvars(phi.im.in); k=k+1)
1896         {
1897            Maps=Maps+list(ratFuncCompratFunc(phi.ratFuncs[j][k],psi.ratFuncs[i]));
1898         }
1899         l=l+1;
1900      }
1901   }
1902
1903   chartmap result=psi.preim,phi.im,Dom,Maps;
1904   kill Dom;
1905   setring R;
1906   return(result);
1907}
1908example
1909{
1910ring R=0,(a,b,c,d),dp;
1911chart U=ideal(ad-bc,ac-bd);
1912Ideal Dom1=ideal(c,d);
1913ratFunc F1=a,c;
1914ratFunc F2=b,c;
1915ratFunc F3=1;
1916ratFunc G1=b,d;
1917ratFunc G2=a,d;
1918ratFunc G3=1;
1919
1920ring S=0,(x,y,z),dp;
1921chart V=xy-z2;
1922chartmap psi=U,V,Dom1,F1,F2,F3,G1,G2,G3;
1923Ideal Dom2=ideal(y,z);
1924ratFunc H1=z,y;
1925ratFunc H2=x,z;
1926
1927ring T=0,(x),dp;
1928chart W=0;
1929chartmap phi=V,W,Dom2,H1,H2;
1930
1931phi*psi;
1932}
1933
1934
1935
1936
1937/////////////////////////////////////////////////////////////////////////////////////
1938//  Morphism Stuff
1939/////////////////////////////////////////////////////////////////////////////////////
1940
1941proc makeMorphism(scheme Preim, scheme Im, list Maps)
1942"USAGE: makeMorphism(Preim,Im,Maps); Preim,Im schemes, Maps list
1943ASSUME: Maps contains lists of the form (phi,i,j), where i,j are integers and phi is a chartmap
1944        from the Preim.cover[i] to Im.cover[j]
1945RETURN: morphism psi from Preim to Im such that psi.chartmaps[Maps[k][2]][Maps[k][3]]=Maps[k][1]
1946        for each k
1947EXAMPLE: example makeMorphism; gives an example"
1948{
1949   morphism result;
1950   def R=basering;
1951
1952   result.preim=Preim;
1953   result.im=Im;
1954
1955   int i,j,k;
1956   for(i=1; i<=size(Preim.cover); i=i+1)
1957   {
1958      result.chartmaps[i]=list();
1959      for(j=1; j<=size(Im.cover); j=j+1)
1960      {
1961         chartmap phi;
1962         phi.preim=Preim.cover[i];
1963         phi.im=Im.cover[j];
1964
1965         def S=Preim.cover[i].in;
1966         setring S;
1967         phi.dom=ideal(0);
1968
1969         result.chartmaps[i][j]=phi;
1970
1971         kill S;
1972         kill phi;
1973      }
1974   }
1975
1976   for(k=1; k<=size(Maps); k=k+1)
1977   {
1978      i=Maps[k][2];
1979      j=Maps[k][3];
1980      result.chartmaps[i][j]=Maps[k][1];
1981   }
1982
1983   setring R;
1984   return(result);
1985}
1986example
1987{"EXAMPLE:"; echo=2;
1988ring R=0,(x,y,z),dp;
1989scheme X=makeProjScheme(xz-y2);
1990ring S=0,(a,b),dp;
1991scheme Y=makeProjScheme(0);
1992
1993def T=X.cover[1].in;
1994setring T;
1995Ideal Dom=1;
1996ratFunc F=y;
1997chartmap phi=X.cover[1],Y.cover[1],Dom,F;
1998list L=list(list(phi,1,1));
1999kill T;
2000
2001def T=X.cover[1].in;
2002setring T;
2003Dom=y;
2004F=1,y;
2005phi=X.cover[1],Y.cover[2],Dom,F;
2006L=L+list(list(phi,1,2));
2007kill T;
2008
2009def T=X.cover[2].in;
2010setring T;
2011Dom=1;
2012F=z;
2013phi=X.cover[2],Y.cover[1],Dom,F;
2014L=L+list(list(phi,2,1));
2015kill T;
2016
2017def T=X.cover[2].in;
2018setring T;
2019Dom=1;
2020F=x;
2021phi=X.cover[2],Y.cover[2],Dom,F;
2022L=L+list(list(phi,2,2));
2023kill T;
2024
2025def T=X.cover[3].in;
2026setring T;
2027Dom=x;
2028F=y,x;
2029phi=X.cover[3],Y.cover[1],Dom,F;
2030L=L+list(list(phi,3,1));
2031kill T;
2032
2033def T=X.cover[3].in;
2034setring T;
2035Dom=1;
2036F=y;
2037phi=X.cover[3],Y.cover[2],Dom,F;
2038L=L+list(list(phi,3,2));
2039
2040
2041morphism psi=X,Y,L;
2042psi.preim.cover;
2043psi.im.cover;
2044psi.chartmaps;
2045}
2046
2047
2048
2049proc MorphismCompMorphism(morphism phi, morphism psi)
2050"USAGE: MorphismCompMorphism(phi,psi); phi,psi morphism
2051RETURN: Composition of phi with psi
2052EXAMPLE: example MorphismCompMorphism; gives an example"
2053{
2054   if(phi.preim!=psi.im)
2055   {
2056      ERROR("Inputs not compatible");
2057   }
2058
2059   def R=basering;
2060   morphism result;
2061
2062   result.preim=psi.preim;
2063   result.im=phi.im;
2064
2065   int i,j,k,l;
2066   for(i=1; i<=size(psi.preim.cover); i=i+1)
2067   {
2068      result.chartmaps[i]=list();
2069      for(j=1; j<=size(phi.im.cover); j=j+1)
2070      {
2071         chartmap f;
2072         f.preim=psi.preim.cover[i];
2073         f.im=phi.im.cover[j];
2074
2075         def S=psi.preim.cover[i].in;
2076         setring S;
2077         f.dom=ideal(0);
2078
2079         for(k=1; k<=size(phi.preim.cover); k=k+1)
2080         {
2081            chartmap g=phi.chartmaps[k][j]*psi.chartmaps[i][k];
2082
2083            if(std(g.dom)!=0)
2084            {
2085               f.dom=f.dom,g.dom;
2086               f.dom=simplify(f.dom,2);
2087               f.ratFuncs=f.ratFuncs+g.ratFuncs;
2088            }
2089
2090            kill g;
2091         }
2092
2093         k=1;
2094         while(k<=size(f.dom))
2095         {
2096             poly p=f.dom[k];
2097             l=1;
2098             while(l<=size(f.dom))
2099             {
2100                if(k==l)
2101                {
2102                   l=l+1;
2103                }else{
2104                if(l<k)
2105                {
2106                   if(reduce(f.dom[l],p)==0)
2107                   {
2108                      f.dom=delete(f.dom,l);
2109                      f.ratFuncs=delete(f.ratFuncs,l);
2110                      k=k-1;
2111                   }else{
2112                      l=l+1;
2113                   }
2114                }else{
2115                   if(reduce(f.dom[l],p)==0)
2116                   {
2117                      f.dom=delete(f.dom,l);
2118                      f.ratFuncs=delete(f.ratFuncs,l);
2119                   }else{
2120                      l=l+1;
2121                   }
2122                }}
2123             }
2124
2125             k=k+1;
2126
2127             kill p;
2128         }
2129
2130         result.chartmaps[i][j]=f;
2131
2132         kill S;
2133         kill f;
2134      }
2135   }
2136
2137   setring R;
2138   return(result);
2139}
2140example
2141{"EXAMPLE:"; echo=2;
2142ring R=0,(x,y,z),dp;
2143scheme X=makeProjScheme(xz-y2);
2144ring S=0,(a,b),dp;
2145scheme Y=makeProjScheme(0);
2146
2147def T=X.cover[1].in;
2148setring T;
2149Ideal Dom=1;
2150ratFunc F=y;
2151chartmap f=X.cover[1],Y.cover[1],Dom,F;
2152list L=list(list(f,1,1));
2153kill T;
2154
2155def T=X.cover[1].in;
2156setring T;
2157Dom=y;
2158F=1,y;
2159f=X.cover[1],Y.cover[2],Dom,F;
2160L=L+list(list(f,1,2));
2161kill T;
2162
2163def T=X.cover[2].in;
2164setring T;
2165Dom=1;
2166F=z;
2167f=X.cover[2],Y.cover[1],Dom,F;
2168L=L+list(list(f,2,1));
2169kill T;
2170
2171def T=X.cover[2].in;
2172setring T;
2173Dom=1;
2174F=x;
2175f=X.cover[2],Y.cover[2],Dom,F;
2176L=L+list(list(f,2,2));
2177kill T;
2178
2179def T=X.cover[3].in;
2180setring T;
2181Dom=x;
2182F=y,x;
2183f=X.cover[3],Y.cover[1],Dom,F;
2184L=L+list(list(f,3,1));
2185kill T;
2186
2187def T=X.cover[3].in;
2188setring T;
2189Dom=1;
2190F=y;
2191f=X.cover[3],Y.cover[2],Dom,F;
2192L=L+list(list(f,3,2));
2193
2194morphism phi=X,Y,L;
2195kill T;
2196
2197
2198def T=Y.cover[1].in;
2199setring T;
2200Dom=1;
2201F=b;
2202ratFunc G=b2;
2203f=Y.cover[1],X.cover[1],Dom,F,G;
2204L=list(list(f,1,1));
2205kill T;
2206
2207def T=Y.cover[1].in;
2208setring T;
2209Dom=b;
2210F=1,b;
2211G=b;
2212f=Y.cover[1],X.cover[2],Dom,F,G;
2213L=L+list(list(f,1,2));
2214kill T;
2215
2216def T=Y.cover[1].in;
2217setring T;
2218Dom=b;
2219F=1,b2;
2220G=1,b;
2221f=Y.cover[1],X.cover[3],Dom,F,G;
2222L=L+list(list(f,1,3));
2223kill T;
2224
2225def T=Y.cover[2].in;
2226setring T;
2227Dom=a;
2228F=1,a;
2229G=1,a2;
2230f=Y.cover[2],X.cover[1],Dom,F,G;
2231L=L+list(list(f,2,1));
2232kill T;
2233
2234def T=Y.cover[2].in;
2235setring T;
2236Dom=a;
2237F=a;
2238G=1,a;
2239f=Y.cover[2],X.cover[2],Dom,F,G;
2240L=L+list(list(f,2,2));
2241kill T;
2242
2243def T=Y.cover[2].in;
2244setring T;
2245Dom=1;
2246F=a2;
2247G=a;
2248f=Y.cover[2],X.cover[3],Dom,F,G;
2249L=L+list(list(f,2,3));
2250
2251morphism psi=Y,X,L;
2252
2253morphism id=psi*phi;
2254id.chartmaps;
2255}
2256
2257
2258proc ChartmapToMorphism(chartmap phi)
2259"USAGE: ChartmapToMorphism(phi); phi chartmap
2260RETURN: morphism corresponding to phi
2261EXAMPLE: example ChartmapToMorphism; gives an example"
2262{
2263   morphism result;
2264
2265   result.preim=ChartToScheme(phi.preim);
2266   result.im=ChartToScheme(phi.im);
2267   result.chartmaps=list(list(phi));
2268
2269   return(result);
2270}
2271example
2272{"EXAMPLE:"; echo=2;
2273ring R=0,(x,y),dp;
2274chart U=0;
2275Ideal Dom=1;
2276list L=list(makeratFunc(x),makeratFunc(y),makeratFunc(1));
2277ring S=0,(a,b,c),dp;
2278chart V=0;
2279chartmap phi=U,V,Dom,L;
2280
2281morphism Phi=ChartmapToMorphism(phi);
2282Phi.preim;
2283Phi.im;
2284Phi.chartmaps;
2285}
2286
2287
2288proc BlowUpMapToMorphism(blowUpMap phi)
2289"USAGE: BlowUpMapToMorphism(phi); phi blowUpMap
2290RETURN: morphism corresponding to phi
2291EXAMPLE: example BlowUpMapToMorphism; gives an example"
2292{
2293   morphism result=phi.blow,phi.base,list();
2294
2295   int i,j;
2296   for(i=1; i<=size(phi.blow.cover); i=i+1)
2297   {
2298      j=phi.maps[i][2];
2299      result.chartmaps[i][j]=phi.maps[i][1];
2300   }
2301
2302   return(result);
2303}
2304example
2305{"EXAMPLE:"; echo=2;
2306ring R=0,(x,y,z),dp;
2307scheme X=makeProjScheme(0);
2308scheme Z=makeProjScheme(ideal(x-z,y-z),X);
2309
2310list L=blowUp(X,Z);
2311
2312morphism Pi=BlowUpMapToMorphism(L[2]);
2313Pi.preim.cover;
2314Pi.im.cover;
2315Pi.chartmaps;
2316}
2317
2318
2319
2320proc InvImage(scheme X, morphism phi)
2321"USAGE: InvImage(x,phi); X scheme, phi morphism
2322RETURN: Inverse image of X under phi
2323EXAMPLE: example InvImage; gives an example"
2324{
2325   def R=basering;
2326   scheme result;
2327
2328   int i,j;
2329   for(i=1; i<=size(phi.preim.cover); i=i+1)
2330   {
2331      def S=phi.preim.cover[i].in;
2332      setring S;
2333      ideal I=1;
2334      Ideal J;
2335      for(j=1; j<=size(phi.im.cover); j=j+1)
2336      {
2337         def T=phi.im.cover[j].in;
2338         setring S;
2339         if(std(phi.chartmaps[i][j].dom)!=0)
2340         {
2341         setring T;
2342         J=ChartmapCompIdeal(phi.chartmaps[i][j],X.cover[j].value);
2343         setring S;
2344         I=intersect(I,J.value);
2345         }
2346
2347         kill T;
2348      }
2349      setring S;
2350      I=std(I+phi.preim.cover[i].value);
2351      chart U=I;
2352      result.cover[i]=U;
2353
2354      kill I;
2355      kill J,S,U;
2356   }
2357
2358   result.maps=phi.preim.maps;
2359
2360   for(i=1; i<=size(result.cover); i=i+1)
2361   {
2362      for(j=1; j<=size(result.cover); j=j+1)
2363      {
2364         result.maps[i][j].preim=result.cover[i];
2365         result.maps[i][j].im=result.cover[j];
2366      }
2367   }
2368
2369   setring R;
2370   return(result);
2371}
2372example
2373{"EXAMPLE:"; echo=2;
2374ring R=0,(x,y,z),dp;
2375scheme X=makeProjScheme(xz-y2);
2376ring S=0,(a,b),dp;
2377scheme Y=makeProjScheme(0);
2378
2379def T=X.cover[1].in;
2380setring T;
2381Ideal Dom=1;
2382ratFunc F=y;
2383chartmap f=X.cover[1],Y.cover[1],Dom,F;
2384list L=list(list(f,1,1));
2385kill T;
2386
2387def T=X.cover[1].in;
2388setring T;
2389Dom=y;
2390F=1,y;
2391f=X.cover[1],Y.cover[2],Dom,F;
2392L=L+list(list(f,1,2));
2393kill T;
2394
2395def T=X.cover[2].in;
2396setring T;
2397Dom=1;
2398F=z;
2399f=X.cover[2],Y.cover[1],Dom,F;
2400L=L+list(list(f,2,1));
2401kill T;
2402
2403def T=X.cover[2].in;
2404setring T;
2405Dom=1;
2406F=x;
2407f=X.cover[2],Y.cover[2],Dom,F;
2408L=L+list(list(f,2,2));
2409kill T;
2410
2411def T=X.cover[3].in;
2412setring T;
2413Dom=x;
2414F=y,x;
2415f=X.cover[3],Y.cover[1],Dom,F;
2416L=L+list(list(f,3,1));
2417kill T;
2418
2419def T=X.cover[3].in;
2420setring T;
2421Dom=1;
2422F=y;
2423f=X.cover[3],Y.cover[2],Dom,F;
2424L=L+list(list(f,3,2));
2425
2426morphism phi=X,Y,L;
2427
2428setring S;
2429scheme Z=makeProjScheme(ideal(3a-2b),Y);
2430
2431scheme Zinv=InvImage(Z,phi);
2432Zinv;
2433}
2434
2435
2436
2437
2438
2439
2440////////////////////////////////////////////////////////////////////////////////////
2441//    Blow-Up
2442////////////////////////////////////////////////////////////////////////////////////
2443
2444proc blowUp(scheme X, scheme C)
2445"USAGE: blowUp(X,C);  X,C scheme
2446ASSUME: C is a subscheme of X
2447RETURN: return a list l=(B,pi,E),
2448          B scheme, is the blow-up of X in the center C
2449          pi blowUpMap, is the map corresponding to the blow-up and contains all the relevant information
2450          E scheme, is the exeptional divisor of the blow-up
2451EXAMPLE: example blowUp;  shows an example
2452"
2453{
2454   scheme B;
2455   blowUpMap pi;
2456   scheme E;
2457   list Bcover=list();
2458   list maps=list();
2459   list eD=list();
2460
2461   pi.base=X;
2462   for(int i=1; i<=size(C.cover); i=i+1)
2463   {
2464       def S=C.cover[i].in;
2465       setring S;
2466       C.cover[i].value=std(C.cover[i].value);
2467       kill S;
2468   }
2469   pi.center=C;
2470
2471   int n=size(X.cover);
2472   for(i=1; i<=n; i=i+1)
2473   {
2474      def R=X.cover[i].in;
2475      setring R;
2476      int m=nvars(R);
2477      C.cover[i].value=simplify(C.cover[i].value,2);
2478      int l=size(C.cover[i].value);
2479      ideal A=C.cover[i].value;
2480      ideal Base=X.cover[i].value;
2481      for(int j=1; j<=l; j=j+1)
2482      {
2483         def S=changevar("x()");
2484
2485         list newvars;
2486         for(int k=1; k<=l; k=k+1)
2487         {
2488            if(k!=j)
2489            {   newvars=newvars+list("y("+string(k)+")");}
2490         }
2491         if(size(newvars)!=0)
2492         {
2493            S=addvarsTo(S,newvars,0);
2494         }
2495
2496         setring S;
2497         chart U=ideal(0);
2498         chartmap phi;
2499         phi.preim=U;
2500         phi.im=X.cover[i];
2501         phi.dom=1;
2502         phi.ratFuncs[1]=list();
2503         for(k=1; k<=m; k=k+1)
2504         {
2505            ratFunc r=var(k);
2506            phi.ratFuncs[1]=phi.ratFuncs[1]+list(r);
2507            kill r;
2508         }
2509
2510
2511
2512         setring S;
2513         ideal I=0;
2514         ideal J=0;
2515         for(k=1; k<=m; k=k+1)
2516         {
2517            I=I+ideal(var(k));
2518         }
2519         map f=R,I;
2520         J=f(A);
2521         I=f(Base);
2522         for(k=1; k<=l; k=k+1)
2523         {
2524            if(k!=j)
2525            {
2526               I=I+ideal(J[j]*y(k)-J[k]);
2527            }
2528         }
2529         U.value=sat(I,ideal(J[j]))[1];
2530         phi.preim=U;
2531         Bcover=Bcover+list(U);
2532         maps=maps+list(list(phi,i,j));
2533         eD=eD+list(makeChart(std(U.value+J)));
2534
2535
2536         kill J,newvars,S,I,U,phi,k;
2537         setring R;
2538      }
2539
2540      kill R,m,l,j;
2541   }
2542
2543   pi.maps=maps;
2544   B.cover=Bcover;
2545   int m=size(B.cover);
2546   int j,r,s,k,l,ll,e,kPreim,kIm,ny,nx,nn;
2547   list L,facts;
2548   ring R;
2549   for(i=1; i<=m; i=i+1)
2550   {
2551      L=list();
2552      for(j=1; j<=m; j=j+1)
2553      {
2554         kPreim=maps[i][3];
2555         kIm=maps[j][3];
2556         r=maps[i][2];
2557         s=maps[j][2];
2558         if(r!=s)
2559         {
2560            chartmap phi;
2561            phi.preim=B.cover[i];
2562            phi.im=B.cover[j];
2563
2564            chartmap psi=X.maps[r][s];
2565            def S=X.cover[r].in;
2566            def SS=X.cover[s].in;
2567            def T=B.cover[i].in;
2568            def TT=B.cover[j].in;
2569
2570            setring T;
2571            ideal I;
2572            for(k=1; k<=nvars(S); k=k+1)
2573            {
2574               I[k]=x(k);
2575            }
2576            map g=S,I;
2577            kill I;
2578
2579            setring S;
2580            poly f,p,q;
2581            for(k=1; k<=size(psi.dom); k=k+1)
2582            {
2583               facts=list();
2584               f=psi.dom[k];
2585               e=sat(C.cover[r].value,ideal(f))[2];
2586               setring SS;
2587               for(l=1; l<=size(C.cover[s].value); l=l+1)
2588               {
2589                  facts[l]=list();
2590                  ratFunc F=PolyCompratFunc(makePoly(C.cover[s].value[l]),psi.ratFuncs[k]);
2591                  setring S;
2592                  list LL=pdivi(F.num*f^e,C.cover[r].value);
2593                  nn=size(C.cover[r].value);
2594                  for(ll=1; ll<=nn; ll=ll+1)
2595                  {
2596                     p=LL[2][ll];
2597                     setring T;
2598                     facts[l][ll]=makePoly(g(p));
2599                     setring S;
2600                  }
2601                  kill LL;
2602                  setring S;
2603                  p=F.den*f^e;
2604                  setring T;
2605                  facts[l][nn+1]=makePoly(g(p));
2606
2607                  kill F;
2608                  setring SS;
2609               }
2610
2611               setring T;
2612               ratFunc F,G;
2613               poly pt;
2614               for(l=1; l<=size(facts[kIm])-1; l=l+1)
2615               {
2616                  if(l!=kPreim)
2617                  {
2618                  pt=pt+facts[kIm][l].value*y(l);
2619                  }else{
2620                  pt=pt+facts[kIm][l].value;
2621                  }
2622               }
2623               pt=pt*g(f);
2624               phi.dom[k]=pt;
2625               phi.ratFuncs[k]=list();
2626               F=pt,facts[kIm][size(facts[kIm])].value;
2627               kill pt;
2628
2629               for(l=1; l<=nvars(SS); l=l+1)
2630               {
2631                  G=ratFuncCompratFunc(psi.ratFuncs[k][l],maps[i][1].ratFuncs[1]);
2632                  phi.ratFuncs[k]=phi.ratFuncs[k]+list(G);
2633               }
2634
2635               for(l=1; l<=size(facts); l=l+1)
2636               {
2637                  if(l!=kIm)
2638                  {
2639                  G=0;
2640                  for(ll=1; ll<=size(facts[l])-1; ll=ll+1)
2641                  {
2642                     if(ll!=kPreim)
2643                     {
2644                     G=G+makeratFunc(facts[l][ll].value*y(ll),facts[l][size(facts[l])].value);
2645                     }else{
2646                     G=G+makeratFunc(facts[l][ll].value,facts[l][size(facts[l])].value);
2647                     }
2648                  }
2649                  phi.ratFuncs[k]=phi.ratFuncs[k]+list(G/F);
2650                  }
2651               }
2652
2653               kill F,G;
2654               setring S;
2655            }
2656
2657            L=L+list(phi);
2658
2659            setring S;
2660            kill f,p,q;
2661
2662            setring T;
2663            kill g;
2664
2665            kill phi,psi,S,SS,T,TT;
2666         }
2667         else
2668         {
2669            def S=X.cover[r].in;
2670            nx=nvars(S);
2671
2672            def T=B.cover[i].in;
2673            chartmap phi;
2674            phi.preim=B.cover[i];
2675            phi.im=B.cover[j];
2676            setring T;
2677            phi.ratFuncs[1]=list();
2678            if(kPreim==kIm)
2679            {
2680               phi.dom=1;
2681               for(k=1;k<=nx;k=k+1)
2682               {
2683                   ratFunc q=x(k);
2684                   phi.ratFuncs[1]=phi.ratFuncs[1]+list(q);
2685                   kill q;
2686               }
2687               ny=nvars(T)-nx+1;
2688               for(k=1;k<=ny;k=k+1)
2689               {
2690                   if(k!=kPreim)
2691                   {
2692                      ratFunc q=y(k);
2693                      phi.ratFuncs[1]=phi.ratFuncs[1]+list(q);
2694                      kill q;
2695                   }
2696               }
2697            }else{
2698               phi.dom=y(kIm);
2699               for(k=1;k<=nx;k=k+1)
2700               {
2701                   ratFunc q=x(k);
2702                   phi.ratFuncs[1]=phi.ratFuncs[1]+list(q);
2703                   kill q;
2704               }
2705               ny=nvars(T)-nx+1;
2706               for(k=1;k<=ny;k=k+1)
2707               {
2708                   if(k!=kIm)
2709                   {
2710                      if(k!=kPreim){
2711                      ratFunc q=y(k),y(kIm);
2712                      }else{
2713                      ratFunc q=1,y(kIm);
2714                      }
2715                      phi.ratFuncs[1]=phi.ratFuncs[1]+list(q);
2716                      kill q;
2717                   }
2718               }
2719            }
2720
2721            L=L+list(phi);
2722            kill S,T,phi;
2723            setring R;
2724         }
2725      }
2726      B.maps=B.maps+list(L);
2727   }
2728
2729   E.cover=eD;
2730   E.maps=B.maps;
2731   for(i=1; i<=size(E.cover); i=i+1)
2732   {
2733      for(j=1; j<=size(E.cover); j=j+1)
2734      {
2735            E.maps[i][j].preim=E.cover[i];
2736            E.maps[i][j].im=E.cover[j];
2737      }
2738   }
2739
2740   pi.blow=B;
2741   pi.exDiv=E;
2742   return(list(B,pi,E));
2743}
2744example
2745{ "EXAMPLE"; echo=2;
2746   ring R=0,(x,y,z),dp;
2747   scheme X=makeProjScheme(ideal(0));
2748   scheme Y=makeProjScheme(ideal(x-z,y-z),X);
2749   list L=blowUp(X,Y);
2750
2751   "Affine cover of the blow-up";
2752   L[1].cover;
2753
2754   "Maps between the affine charts";
2755   L[1].maps;
2756
2757   "Exeptional divisor";
2758   L[3].cover;
2759}
2760
2761
2762
2763proc strictTransform(Y,blowUpMap pi)
2764"USAGE: strinctTransform(Y,pi);
2765RETURN: the strict transform of subscheme Y under the blow-up map pi
2766EXAMPLE: example strictTransform, shows an example"
2767{
2768   if(typeof(Y)=="scheme")
2769   {
2770      def B=basering;
2771      scheme result;
2772      result.maps=pi.blow.maps;
2773      int n=size(pi.blow.cover);
2774      for(int k=1; k<=n; k=k+1)
2775      {
2776         int m=pi.maps[k][2];
2777         def R=Y.cover[m].in;
2778         def S=pi.blow.cover[k].in;
2779
2780         setring R;
2781         ideal base=Y.cover[m].value;
2782         int nx=nvars(R);
2783
2784         setring S;
2785         ideal I=var(1);
2786         for(int i=2; i<=nx; i=i+1)
2787         {
2788            I=I,ideal(var(i));
2789         }
2790         map phi=R,I;
2791         ideal J=phi(base);
2792         J=J+pi.blow.cover[k].value;
2793
2794         ideal E=pi.exDiv.cover[k].value;
2795
2796         J=sat(J,E)[1];
2797         chart U=J;
2798         result.cover=result.cover+list(U);
2799
2800         kill S,i,I,E,nx,U,J,phi,m;
2801         setring R;
2802         kill base;
2803         kill R;
2804      }
2805
2806      for(int i=1; i<=n; i=i+1)
2807      {
2808         for(int j=1; j<=n; j=j+1)
2809         {
2810            result.maps[i][j].preim=result.cover[i];
2811            result.maps[i][j].im=result.cover[j];
2812         }
2813         kill j;
2814      }
2815
2816      setring B;
2817      return(result);
2818   }
2819}
2820example
2821{
2822   "EXAMPLE:"; echo=2;
2823   ring R=0,(x,y),dp;
2824   scheme X=makeAffineScheme(ideal(0));
2825   scheme Y=makeAffineScheme(ideal(x2-y3-y2));
2826   scheme Z=makeAffineScheme(ideal(x,y));
2827   list L=blowUp(X,Z);
2828   scheme blowY=strictTransform(Y,L[2]);
2829   blowY.maps;
2830}
2831
2832
2833
2834
2835
2836
2837
2838
2839
2840
2841
2842
2843
2844
2845
2846//////////////////////////////////////////////////////////////////////////////////
2847//  Smoothness Test
2848//////////////////////////////////////////////////////////////////////////////////
2849
2850
2851// Used in CompleteIntersectionCover and HybrSmoothTest
2852static proc findMaxMinors(ideal I, matrix J)
2853"USAGE: findMaxMinors(I,J); I ideal, J matrix
2854ASSUME: J is Jacobi matrix of II
2855RETURN: A list of lists of the form (h,A,iv) where h is a maximal minor of J and iv contains
2856        the colomns used in the formation of h and A is the matrix of cofactors of the matrix
2857        forming h
2858"
2859{
2860   list Minors;
2861   ideal Inew=std(I);
2862   poly h;
2863   int n=ncols(J);
2864   int r=nrows(J);
2865   if(r>1)
2866   {
2867   matrix M[r][r];
2868   matrix A[r][r];
2869   matrix B[r][r];
2870   intvec iv=1..r;
2871   int done;
2872   int k,l,c,i,j;
2873   k=r;
2874   c=0;
2875   while((k>0)&&(Inew!=1))
2876   {
2877      M=submat(J,1..r,iv);
2878      h=det(M);
2879      if(h!=0)
2880      {
2881         c=c+1;
2882         for(i=1; i<=r; i=i+1)
2883         {
2884         for(j=1; j<=r; j=j+1)
2885         {
2886            B=permrow(M,i,r);
2887            B=permcol(B,j,r);
2888            A[i,j]=det(submat(B,1..r-1,1..r-1))*(-1)^(i+j);
2889         }}
2890         A=transpose(A);
2891         Minors[c]=list(h,A,iv);
2892         Inew=std(radical(Inew+h));
2893      }
2894
2895      k=r-1;
2896      if(iv[r]<n)
2897      {
2898         iv[r]=iv[r]+1;
2899      }else{
2900         done=0;
2901         while(!done)
2902         {
2903            if(iv[k]>=(iv[k+1]-1))
2904            {
2905               k=k-1;
2906               if(k<=0)
2907               { done=1;}
2908            }else{ done=1;}
2909         }
2910         if(k>0)
2911         {
2912            iv[k]=iv[k]+1;
2913            l=1;
2914            while(k+l<=r)
2915            {
2916               iv[k+l]=iv[k+l-1]+1;
2917               l=l+1;
2918            }
2919         }
2920      }
2921   }
2922   return(Minors);
2923   }else{
2924   matrix A[r][r];
2925   int i,c;
2926   intvec iv;
2927   c=0;
2928   i=1;
2929   while((i<=n)&&(Inew!=1))
2930   {
2931      iv=i;
2932      h=J[1,i];
2933      if(h!=0)
2934      {
2935         c=c+1;
2936         A[1,1]=1;
2937         Minors[c]=list(h,A,iv);
2938         Inew=std(radical(Inew+h));
2939      }
2940      i=i+1;
2941   }
2942   return(Minors);
2943   }
2944}
2945
2946
2947proc CompleteIntersectionCover(list #)
2948"USAGE: CompleteIntersectionCover(X,W); X,W schemes
2949ASSUME: X is a subscheme of W, X and W are equidimensional and smooth and the charts of W are
2950        complete intersections
2951RETURN: A list L where:
2952            L[1] is a scheme isomorphic to X where each chart is a complete intersection
2953            L[2] is the morphism representing the isomorphism X->L[1]
2954            L[3] is the morphism representing the isomorphism L[1]->X
2955EXAMPLE: example CompleteIntersectionCover; gives an example"
2956{
2957if((size(#)==1)&&(typeof(#[1])=="ideal"))
2958{
2959   def R=basering;
2960   ideal I=#[1];
2961   scheme resultScheme;
2962   morphism resultMap,resultInvMap;
2963   if((std(I)==0)||(std(I)==1))
2964   {
2965      resultScheme=makeAffineScheme(std(I));
2966      resultMap.preim=resultScheme;
2967      resultMap.im=resultScheme;
2968      resultMap.chartmaps=resultScheme.maps;
2969      resultInvMap.preim=resultScheme;
2970      resultInvMap.im=resultScheme;
2971      resultInvMap.chartmaps=resultScheme.maps;
2972      return(list(resultScheme,resultMap,resultInvMap));
2973   }
2974   else
2975   {
2976      matrix J=jacob(I);
2977      int i,j,n,k,l,ci,cj;
2978      scheme X,W;
2979      list L,K;
2980      ideal H=std(radical(I));
2981      n=0;
2982      for(i=1; i<=size(I); i=i+1)
2983      {
2984         for(j=1; j<=nvars(R); j=j+1)
2985         {
2986            if(H!=1)
2987            {
2988            if(reduce(J[i,j],H)!=0)
2989            {
2990            n=n+1;
2991            L[n]=J[i,j];
2992            K[n]=i;
2993            H=std(radical(H+J[i,j]));
2994            }}
2995         }
2996      }
2997
2998      list CoverSchemes,PreCoverCharts,PreCoverMaps,CoverMaps,CoverInvMaps;
2999      for(k=1; k<=n; k=k+1)
3000      {
3001         poly h=L[k];
3002         def S=changevar("x()",R);
3003         S=extendring(1,"t","dp",0,S);
3004         setring S;
3005         ideal VarR=var(2);
3006         for(j=2; j<=nvars(R); j=j+1)
3007         {
3008            VarR=VarR,var(j+1);
3009         }
3010         map phi=R,VarR;
3011         i=K[k];
3012         poly hnew=phi(h);
3013         poly g=t*hnew-1;
3014
3015         ideal Inew=phi(I);
3016         poly f=Inew[i];
3017         Inew[i]=Inew[1];
3018         Inew[1]=f;
3019         Inew=g,Inew;
3020
3021         W=makeAffineScheme(ideal(g,f));
3022         X=makeAffineScheme(Inew,W);
3023
3024         PreCoverCharts[k]=X.cover[1];
3025         list LL=CompleteIntersectionCover(Inew,ideal(g,f));
3026         CoverSchemes[k]=LL[1];
3027         CoverMaps[k]=LL[2];
3028         CoverInvMaps[k]=LL[3];
3029         kill LL;
3030
3031         kill phi,hnew,Inew,g,VarR;
3032         kill S;
3033         setring R;
3034         kill h;
3035      }
3036
3037      for(i=1; i<=n; i=i+1)
3038      {
3039         list Maps;
3040         poly h;
3041         chart Preim=PreCoverCharts[i];
3042         def S=Preim.in;
3043         setring S;
3044         ideal VarR=var(2);
3045         for(j=2; j<=nvars(R); j=j+1)
3046         {
3047            VarR=VarR,var(j+1);
3048         }
3049         map psi=R,VarR;
3050         setring R;
3051
3052         for(j=1; j<=n; j=j+1)
3053         {
3054            chartmap phi;
3055            phi.preim=Preim;
3056            phi.im=PreCoverCharts[j];
3057            h=L[j];
3058            setring S;
3059            poly hnew=psi(h);
3060            phi.dom=hnew;
3061            list RatFuncs;
3062            RatFuncs[1]=makeratFunc(1,hnew);
3063            for(k=1; k<=nvars(R); k=k+1)
3064            {
3065               RatFuncs[k+1]=makeratFunc(var(k+1));
3066            }
3067            phi.ratFuncs[1]=RatFuncs;
3068            Maps[j]=phi;
3069
3070            kill hnew,RatFuncs;
3071            setring R;
3072            kill phi;
3073         }
3074         PreCoverMaps[i]=Maps;
3075
3076         setring S;
3077         kill psi;
3078         kill Maps,Preim,S;
3079         setring R;
3080         kill h;
3081      }
3082
3083      X=makeAffineScheme(I);
3084      list Maps;
3085      resultMap.preim=X;
3086      resultInvMap.im=X;
3087      for(i=1; i<=n; i=i+1)
3088      {
3089         resultScheme.cover=resultScheme.cover+CoverSchemes[i].cover;
3090
3091         chartmap phi,psi;
3092         phi.preim=X.cover[1]; phi.im=PreCoverCharts[i];
3093         psi.im=X.cover[1]; psi.preim=PreCoverCharts[i];
3094         list phiRatFuncs;
3095         list psiRatFuncs;
3096
3097         phiRatFuncs[1]=makeratFunc(1,L[i]);
3098         for(j=1; j<=nvars(R); j=j+1)
3099         {
3100            phiRatFuncs=phiRatFuncs+list(makeratFunc(var(j)));
3101         }
3102         phi.dom=ideal(L[i]);
3103         phi.ratFuncs[1]=phiRatFuncs;
3104
3105         def S=PreCoverCharts[i].in;
3106         setring S;
3107         for(j=1; j<=nvars(R); j=j+1)
3108         {
3109            psiRatFuncs=psiRatFuncs+list(makeratFunc(var(j+1)));
3110         }
3111         psi.dom=1;
3112         psi.ratFuncs[1]=psiRatFuncs;
3113         setring R;
3114
3115
3116         for(k=1; k<=size(CoverSchemes[i].cover); k=k+1)
3117         {
3118            Maps=Maps+list(CoverMaps[i].chartmaps[1][k]*phi);
3119            resultInvMap.chartmaps=resultInvMap.chartmaps+list(list(psi*CoverInvMaps[i].chartmaps[k][1]));
3120         }
3121         kill phi,psi,phiRatFuncs,psiRatFuncs,S;
3122      }
3123      resultMap.chartmaps=list(Maps);
3124
3125
3126      ci=0;
3127      chartmap phi;
3128      for(i=1; i<=n; i=i+1)
3129      {
3130      for(k=1; k<=size(CoverSchemes[i].cover); k=k+1)
3131      {
3132         ci=ci+1;
3133         cj=0;
3134         Maps=list();
3135         for(j=1; j<=n; j=j+1)
3136         {
3137         for(l=1; l<=size(CoverSchemes[j].cover); l=l+1)
3138         {
3139            cj=cj+1;
3140            if(i==j)
3141            {
3142               Maps[cj]=CoverSchemes[i].maps[k][l];
3143            }else{
3144               phi=CoverMaps[j].chartmaps[1][l];
3145               phi=phi*PreCoverMaps[i][j];
3146               phi=phi*CoverInvMaps[i].chartmaps[k][1];
3147               Maps[cj]=phi;
3148            }
3149         }}
3150         resultScheme.maps[ci]=Maps;
3151      }}
3152      resultMap.im=resultScheme;
3153      resultInvMap.preim=resultScheme;
3154
3155      return(list(resultScheme,resultMap,resultInvMap));
3156   }
3157}
3158
3159if((size(#)==2)&&(typeof(#[1])=="ideal")&&(typeof(#[2])=="ideal"))
3160{
3161   def R=basering;
3162   ideal IX=#[1];
3163   ideal IW=#[2];
3164   if(std(IW)==0)
3165   {   return(CompleteIntersectionCover(IX));}
3166   int k;
3167   int r=size(IW);
3168   int s=size(IX);
3169   if(r>s)
3170   {   ERROR("Incompatible Inputs");}
3171   for(k=1; k<=r; k=k+1)
3172   {   if(IW[k]!=IX[k]){ERROR("Incompatible Inputs");}}
3173   scheme resultScheme;
3174   morphism resultMap,resultInvMap;
3175   if(dim(std(IX))==dim(std(IW)))
3176   {
3177      resultScheme=makeAffineScheme(IW);
3178      resultMap.preim=resultScheme;
3179      resultMap.im=resultScheme;
3180      resultMap.chartmaps=resultScheme.maps;
3181      resultInvMap.preim=resultScheme;
3182      resultInvMap.im=resultScheme;
3183      resultInvMap.chartmaps=resultScheme.maps;
3184      return(list(resultScheme,resultMap,resultInvMap));
3185   }
3186   else
3187   {
3188      matrix JX=jacob(IX);
3189      matrix JW=jacob(IW);
3190      int i,j,n,m,l,ci,cj;
3191      scheme X,W;
3192      list L,K;
3193
3194      list Minors=findMaxMinors(IW,JW);
3195      m=size(Minors);
3196      n=0;
3197      for(k=1; k<=m; k=k+1)
3198      {
3199         matrix v[s][1];
3200         for(i=1; i<=s; i=i+1)
3201         {
3202            v[i,1]=IX[i];
3203         }
3204         matrix M=dsum(Minors[k][2],diag(Minors[k][1],s-r));
3205         ideal I=M*v;
3206         matrix J=jacob(I);
3207         intvec iv=Minors[k][3];
3208         for(i=r+1; i<=s; i=i+1)
3209         {
3210            for(j=1; j<=r; j=j+1)
3211            {
3212               poly g=-diff(IX[i],var(iv[j]));
3213               J=addrow(J,j,g,i);
3214               kill g;
3215            }
3216         }
3217
3218         ideal H=std(radical(IW));
3219         for(i=r+1; i<=s; i=i+1)
3220         {
3221         l=1;
3222         for(j=1; j<=nvars(R); j=j+1)
3223         {
3224            if(reduce(Minors[k][1],H)!=0)
3225            {
3226            if(j!=iv[l])
3227            {
3228            if(reduce(J[i,j],H)!=0)
3229            {
3230            n=n+1;
3231            L[n]=J[i,j]*Minors[k][1];
3232            K[n]=i;
3233            H=std(radical(H+J[i,j]));
3234            }}
3235            else{
3236            if(l<r){l=l+1;}
3237            }}
3238         }
3239         }
3240
3241         kill v,M,I,J,iv,H;
3242      }
3243
3244      list CoverSchemes,PreCoverCharts,PreCoverMaps,CoverMaps,CoverInvMaps;
3245      for(k=1; k<=n; k=k+1)
3246      {
3247         poly h=L[k];
3248         def S=changevar("x()",R);
3249         S=extendring(1,"t","dp",0,S);
3250         setring S;
3251         ideal VarR=var(2);
3252         for(j=2; j<=nvars(R); j=j+1)
3253         {
3254            VarR=VarR,var(j+1);
3255         }
3256         map phi=R,VarR;
3257         i=K[k];
3258         poly hnew=phi(h);
3259         poly g=t*hnew-1;
3260         poly f;
3261
3262         ideal IXnew=phi(IX);
3263         ideal IWnew=phi(IW);
3264         if(std(IWnew)!=0)
3265         {
3266         f=IXnew[i];
3267         IXnew[i]=IXnew[r+1];
3268         IXnew[r+1]=f;
3269         IXnew=g,IXnew;
3270         IWnew=g,IWnew,f;
3271         }else{
3272         f=IXnew[i];
3273         IXnew[i]=IXnew[1];
3274         IXnew[1]=f;
3275         IXnew=g,IXnew;
3276         IWnew=g,f;
3277         }
3278
3279         W=makeAffineScheme(IWnew);
3280         X=makeAffineScheme(IXnew,W);
3281
3282         PreCoverCharts[k]=X.cover[1];
3283         list LL=CompleteIntersectionCover(IXnew,IWnew);
3284         CoverSchemes[k]=LL[1];
3285         CoverMaps[k]=LL[2];
3286         CoverInvMaps[k]=LL[3];
3287         kill LL;
3288
3289         kill phi,hnew,VarR,f,g,IXnew,IWnew;
3290         setring R;
3291         kill S;
3292         kill h;
3293      }
3294
3295      for(i=1; i<=n; i=i+1)
3296      {
3297         list Maps;
3298         poly h;
3299         chart Preim=PreCoverCharts[i];
3300         def S=Preim.in;
3301         setring S;
3302         ideal VarR=var(2);
3303         for(j=2; j<=nvars(R); j=j+1)
3304         {
3305            VarR=VarR,var(j+1);
3306         }
3307         map psi=R,VarR;
3308         setring R;
3309
3310         for(j=1; j<=n; j=j+1)
3311         {
3312            chartmap phi;
3313            phi.preim=Preim;
3314            phi.im=PreCoverCharts[j];
3315            h=L[j];
3316            setring S;
3317            poly hnew=psi(h);
3318            phi.dom=hnew;
3319            list RatFuncs;
3320            RatFuncs[1]=makeratFunc(1,hnew);
3321            for(k=1; k<=nvars(R); k=k+1)
3322            {
3323               RatFuncs[k+1]=makeratFunc(var(k+1));
3324            }
3325            phi.ratFuncs[1]=RatFuncs;
3326            Maps[j]=phi;
3327
3328            kill hnew,RatFuncs;
3329            setring R;
3330            kill phi;
3331         }
3332         PreCoverMaps[i]=Maps;
3333         setring S;
3334         kill VarR,psi;
3335         setring R;
3336         kill Maps,Preim,S,h;
3337      }
3338
3339      X=makeAffineScheme(IX);
3340      list Maps;
3341      resultMap.preim=X;
3342      resultInvMap.im=X;
3343      for(i=1; i<=n; i=i+1)
3344      {
3345         resultScheme.cover=resultScheme.cover+CoverSchemes[i].cover;
3346
3347         chartmap phi,psi;
3348         phi.preim=X.cover[1]; phi.im=PreCoverCharts[i];
3349         psi.im=X.cover[1]; psi.preim=PreCoverCharts[i];
3350         list phiRatFuncs;
3351         list psiRatFuncs;
3352
3353         phiRatFuncs[1]=makeratFunc(1,L[i]);
3354         for(j=1; j<=nvars(R); j=j+1)
3355         {
3356            phiRatFuncs=phiRatFuncs+list(makeratFunc(var(j)));
3357         }
3358         phi.dom=ideal(L[i]);
3359         phi.ratFuncs[1]=phiRatFuncs;
3360
3361         def S=PreCoverCharts[i].in;
3362         setring S;
3363         for(j=1; j<=nvars(R); j=j+1)
3364         {
3365            psiRatFuncs=psiRatFuncs+list(makeratFunc(var(j+1)));
3366         }
3367         psi.dom=1;
3368         psi.ratFuncs[1]=psiRatFuncs;
3369         setring R;
3370
3371         for(k=1; k<=size(CoverSchemes[i].cover); k=k+1)
3372         {
3373            Maps=Maps+list(CoverMaps[i].chartmaps[1][k]*phi);
3374            resultInvMap.chartmaps=resultInvMap.chartmaps+list(list(psi*CoverInvMaps[i].chartmaps[k][1]));
3375         }
3376         kill phi,psi,phiRatFuncs,psiRatFuncs,S;
3377      }
3378      resultMap.chartmaps=list(Maps);
3379
3380      ci=0;
3381      chartmap phi;
3382      for(i=1; i<=n; i=i+1)
3383      {
3384      for(k=1; k<=size(CoverSchemes[i].cover); k=k+1)
3385      {
3386         ci=ci+1;
3387         cj=0;
3388         Maps=list();
3389         for(j=1; j<=n; j=j+1)
3390         {
3391         for(l=1; l<=size(CoverSchemes[j].cover); l=l+1)
3392         {
3393            cj=cj+1;
3394            if(i==j)
3395            {
3396               Maps[cj]=CoverSchemes[i].maps[k][l];
3397            }else{
3398               phi=CoverMaps[j].chartmaps[1][l];
3399               phi=phi*PreCoverMaps[i][j];
3400               phi=phi*CoverInvMaps[i].chartmaps[k][1];
3401               Maps[cj]=phi;
3402            }
3403         }}
3404         resultScheme.maps[ci]=Maps;
3405      }}
3406      resultMap.im=resultScheme;
3407      resultInvMap.preim=resultScheme;
3408
3409      return(list(resultScheme,resultMap,resultInvMap));
3410   }
3411}
3412
3413if((size(#)==2)&&(typeof(#[1])=="scheme")&&(typeof(#[2])=="scheme"))
3414{
3415   def R=basering;
3416   scheme X=#[1];
3417   scheme W=#[2];
3418   scheme resultScheme;
3419   int i,j,k,l;
3420   int n=size(X.cover);
3421   morphism resultMap,resultInvMap;
3422
3423   list CoverSchemes,CoverMaps,CoverInvMaps;
3424   for(i=1; i<=n; i=i+1)
3425   {
3426      def S=X.cover[i].in;
3427      setring S;
3428      ideal IX=X.cover[i].value;
3429      ideal IW=W.cover[i].value;
3430      list LL=CompleteIntersectionCover(IX,IW);
3431      CoverSchemes[i]=LL[1];
3432      CoverMaps[i]=LL[2];
3433      CoverInvMaps[i]=LL[3];
3434
3435      kill IX,IW;
3436      setring R;
3437      kill S,LL;
3438   }
3439
3440
3441   list Maps;
3442   resultMap.preim=X;
3443   resultInvMap.im=X;
3444   for(i=1; i<=n; i=i+1)
3445   {
3446      resultScheme.cover=resultScheme.cover+CoverSchemes[i].cover;
3447   }
3448
3449   for(i=1; i<=n; i=i+1)
3450   {
3451      Maps=list();
3452      for(j=1; j<=n; j=j+1)
3453      {
3454      for(k=1; k<=size(CoverSchemes[j].cover); k=k+1)
3455      {
3456         if(i!=j)
3457         {
3458            chartmap phi;
3459            phi.preim=X.cover[i];
3460            phi.im=CoverSchemes[j].cover[k];
3461
3462            def S=X.cover[i].in;
3463            setring S;
3464            phi.dom=0;
3465
3466            def SS=CoverSchemes[j].cover[k].in;
3467            phi.ratFuncs[1]=list();
3468            for(l=1; l<=nvars(SS); l=l+1)
3469            {
3470               phi.ratFuncs[1]=phi.ratFuncs[1]+list(makeratFunc(0));
3471            }
3472            Maps=Maps+list(phi);
3473            kill phi,S,SS;
3474            setring R;
3475         }else{
3476            Maps=Maps+list(CoverMaps[i].chartmaps[1][k]);
3477         }
3478      }}
3479      resultMap.chartmaps[i]=Maps;
3480   }
3481
3482   for(j=1; j<=n; j=j+1)
3483   {
3484   for(k=1; k<=size(CoverSchemes[j].cover); k=k+1)
3485   {
3486      Maps=list();
3487      for(i=1; i<=n; i=i+1)
3488      {
3489         if(i!=j)
3490         {
3491            chartmap phi;
3492            phi.im=X.cover[i];
3493            phi.preim=CoverSchemes[j].cover[k];
3494
3495            def S=CoverSchemes[j].cover[k].in;
3496            setring S;
3497            phi.dom=0;
3498
3499            def SS=X.cover[i].in;
3500            phi.ratFuncs[1]=list();
3501            for(l=1; l<=nvars(SS); l=l+1)
3502            {
3503               phi.ratFuncs[1]=phi.ratFuncs[1]+list(makeratFunc(0));
3504            }
3505            Maps=Maps+list(phi);
3506            kill phi,S,SS;
3507            setring R;
3508         }else{
3509            Maps=Maps+list(CoverInvMaps[i].chartmaps[k][1]);
3510         }
3511      }
3512      resultInvMap.chartmaps=resultInvMap.chartmaps+list(Maps);
3513   }}
3514
3515   for(i=1; i<=n; i=i+1)
3516   {
3517   for(k=1; k<=size(CoverSchemes[i].cover); k=k+1)
3518   {
3519      Maps=list();
3520      for(j=1; j<=n; j=j+1)
3521      {
3522      for(l=1; l<=size(CoverSchemes[j].cover); l=l+1)
3523      {
3524         chartmap phi=CoverMaps[j].chartmaps[1][l];
3525         phi=phi*X.maps[i][j];
3526         phi=phi*CoverInvMaps[i].chartmaps[k][1];
3527         Maps=Maps+list(phi);
3528         kill phi;
3529      }}
3530      resultScheme.maps=resultScheme.maps+list(Maps);
3531   }}
3532   resultMap.im=resultScheme;
3533   resultInvMap.preim=resultScheme;
3534
3535   return(list(resultScheme,resultMap,resultInvMap));
3536}
3537}
3538example
3539{"EXAMPLE:"; echo=2;
3540ring R=0,(x,y,z),dp;
3541scheme X=makeProjScheme(0);
3542scheme Z=makeProjScheme(ideal(x-z,y-z)*ideal(x+z,y-z),X);
3543
3544list LL=CompleteIntersectionCover(Z,X);
3545LL[1];
3546}
3547
3548
3549
3550proc HybrSmoothTest(list #)
3551"USAGE: HybrSmoothTest(X[,W]); X,W schemes
3552RETURN: 1 if X is smooth, 0 otherwise
3553EXAMPLE: example HybrSmoothTest; shows an example"
3554{
3555if((size(#)==2)&&(typeof(#[1])=="chart")&&(typeof(#[2])=="chart"))
3556{
3557   int result=1;
3558   def T=basering;
3559   if(#[1].in!=#[2].in)
3560   {   ERROR("Incompatible Inputs");}
3561   def R=#[1].in;
3562   setring R;
3563   ideal IX=#[1].value;
3564   ideal IW=#[2].value;
3565   int k;
3566   int r=size(IW);
3567   int s=size(IX);
3568   if(r>s)
3569   {   ERROR("Incompatible Inputs");}
3570   for(k=1; k<=r; k=k+1)
3571   {   if(IW[k]!=IX[k]){~:ERROR("Incompatible Inputs");}}
3572   if(dim(std(IX))==dim(std(IW)))
3573   {
3574      if((std(reduce(std(IX),std(IW)))==0)&&(std(reduce(std(IW),std(IX)))==0))
3575      {   result=1;}
3576      else
3577      {   result=0;}
3578      setring T;
3579      return(result);
3580   }
3581   else
3582   {
3583      matrix JX=jacob(IX);
3584      matrix JW=jacob(IW);
3585      int i,j,n,m,l,ci,cj;
3586      list L,K;
3587      if(std(IW)!=0)
3588      {
3589      list Minors=findMaxMinors(IW,JW);
3590      m=size(Minors);
3591      n=0;
3592      for(k=1; k<=m; k=k+1)
3593      {
3594         matrix v[s][1];
3595         for(i=1; i<=s; i=i+1)
3596         {
3597            v[i,1]=IX[i];
3598         }
3599         matrix M=dsum(Minors[k][2],diag(Minors[k][1],s-r));
3600         ideal I=M*v;
3601         matrix J=jacob(I);
3602         intvec iv=Minors[k][3];
3603         for(i=r+1; i<=s; i=i+1)
3604         {
3605            for(j=1; j<=r; j=j+1)
3606            {
3607               poly g=-diff(IX[i],var(iv[j]));
3608               J=addrow(J,j,g,i);
3609               kill g;
3610            }
3611         }
3612
3613         ideal H=std(radical(IW));
3614         for(i=r+1; i<=s; i=i+1)
3615         {
3616         l=1;
3617         for(j=1; j<=nvars(R); j=j+1)
3618         {
3619            if(reduce(Minors[k][1],H)!=0)
3620            {
3621            if(j!=iv[l])
3622            {
3623            if(reduce(J[i,j],H)!=0)
3624            {
3625            n=n+1;
3626            L[n]=J[i,j]*Minors[k][1];
3627            K[n]=i;
3628            H=std(radical(H+J[i,j]));
3629            }}
3630            else{
3631            if(l<r){l=l+1;}
3632            }}
3633         }
3634         }
3635
3636         if(std(H)!=1)
3637         {   result=0;
3638             return(result);
3639         }
3640
3641         kill v,M,I,J,iv,H;
3642      }
3643      }else{
3644      ideal H=std(radical(IX));
3645      matrix J=jacob(IX);
3646      n=0;
3647      for(i=1; i<=size(IX); i=i+1)
3648      {
3649         for(j=1; j<=nvars(R); j=j+1)
3650         {
3651            if(H!=1)
3652            {
3653            if(reduce(J[i,j],H)!=0)
3654            {
3655            n=n+1;
3656            L[n]=J[i,j];
3657            K[n]=i;
3658            H=std(radical(H+J[i,j]));
3659            }}
3660         }
3661      }
3662      if(H!=1)
3663      {   result=0;
3664          return(result);
3665      }
3666      }
3667
3668
3669      for(k=1; k<=n; k=k+1)
3670      {
3671         poly h=L[k];
3672         def S=changevar("x()",R);
3673         S=extendring(1,"t","dp",0,S);
3674         setring S;
3675         setring S;
3676         ideal VarR=var(2);
3677         for(j=2; j<=nvars(R); j=j+1)
3678         {
3679            VarR=VarR,var(j+1);
3680         }
3681         map phi=R,VarR;
3682         i=K[k];
3683         poly hnew=phi(h);
3684         poly g=t*hnew-1;
3685         poly f;
3686
3687         ideal IXnew=phi(IX);
3688         ideal IWnew=phi(IW);
3689         if(std(IWnew)!=0)
3690         {
3691         f=IXnew[i];
3692         IXnew[i]=IXnew[r+1];
3693         IXnew[r+1]=f;
3694         IXnew=g,IXnew;
3695         IWnew=g,IWnew,f;
3696         }else{
3697         f=IXnew[i];
3698         IXnew[i]=IXnew[1];
3699         IXnew[1]=f;
3700         IXnew=g,IXnew;
3701         IWnew=g,f;
3702         }
3703
3704         chart X=IXnew;
3705         chart W=IWnew;
3706         if(!HybrSmoothTest(X,W))
3707         {   result=0;
3708             return(result);
3709         }
3710
3711
3712         kill phi,hnew,VarR,f,g,IXnew,IWnew,X,W;
3713         kill S;
3714         setring R;
3715         kill h;
3716       }
3717
3718       setring T;
3719       return(result);
3720    }
3721}
3722
3723if((size(#)==1)&&(typeof(#[1])=="scheme"))
3724{
3725   def R=basering;
3726   int result=1;
3727   scheme X=#[1];
3728   for(int n=1; n<=size(X.cover); n=n+1)
3729   {
3730      def S=X.cover[n].in;
3731      setring S;
3732      chart W=0;
3733      setring R;
3734      kill S;
3735      if(!HybrSmoothTest(X.cover[n],W))
3736      {   result=0;}
3737      kill W;
3738   }
3739   return(result);
3740}
3741
3742if((size(#)==2)&&(typeof(#[1])=="scheme")&&(typeof(#[2])=="scheme"))
3743{
3744   int result=1;
3745   scheme X=#[1];
3746   scheme W=#[2];
3747   if(size(W.cover)!=size(X.cover))
3748   {   ERROR("Incompatible Inputs");}
3749   for(int n=1; n<=size(X.cover); n=n+1)
3750   {
3751      if(!HybrSmoothTest(X.cover[n],W.cover[n]))
3752      {   result=0;}
3753   }
3754   return(result);
3755}
3756}
3757example
3758{"EXAMPLE:"; echo=2;
3759
3760ring R=0,(x,y),dp;
3761"Example of a smooth scheme";
3762scheme X=makeAffineScheme(ideal(y-x2));
3763HybrSmoothTest(X);
3764
3765"Example of a scheme that isn't smooth";
3766scheme Y=makeAffineScheme(ideal(y2-x3));
3767HybrSmoothTest(Y);
3768
3769"Blowing up resolves singularities";
3770scheme C=makeAffineScheme(ideal(x,y),Y);
3771scheme Z=blowUp(Y,C)[1];
3772HybrSmoothTest(Z);
3773}
3774
3775
3776
3777proc isSmooth(scheme X)
3778"USAGE: isSmooth(X); scheme X
3779RETURN: 1 if X is smooth, 0 otherwise
3780EXAMPLE: example isSmooth; gives an example"
3781{
3782   def R=basering;
3783   int i,j,k;
3784
3785   if(!isReduced(X))
3786   {   return(0);}
3787
3788   list Dec=IrredDec(X);
3789   for(i=1; i<=size(X.cover); i=i+1)
3790   {
3791      def S=X.cover[i].in;
3792      setring S;
3793      for(j=1; j<=size(Dec)-1; j=j+1)
3794      {
3795      for(k=j+1; k<=size(Dec); k=k+1)
3796      {
3797         if(std(Dec[j][1].cover[i].value+Dec[k][1].cover[i].value)!=1)
3798         {
3799            return(0);
3800         }
3801      }}
3802      kill S;
3803   }
3804
3805   setring R;
3806   for(j=1; j<=size(Dec); j=j+1)
3807   {
3808      if(!HybrSmoothTest(Dec[j][1]))
3809      {
3810         return(0);
3811      }
3812   }
3813
3814   setring R;
3815   return(1);
3816}
3817example
3818{"EXAMPLE:"; echo=2;
3819ring R=0,(x,y,z),dp;
3820scheme X=makeProjScheme(x2y3z5);
3821isSmooth(X);
3822
3823X=makeProjScheme(ideal(xyz-x3));
3824isSmooth(X);
3825
3826X=makeAffineScheme(ideal(x2y-x));
3827isSmooth(X);
3828
3829X=makeProjScheme(y2x-x3-x2z);
3830isSmooth(X);
3831
3832scheme Z=makeProjScheme(ideal(x,y),X);
3833list L=blowUp(X,Z);
3834isSmooth(L[1]);
3835}
Note: See TracBrowser for help on using the repository browser.