Changeset 555878 in git
 Timestamp:
 Feb 5, 2010, 1:14:43 PM (13 years ago)
 Branches:
 (u'spielwiese', '828514cf6e480e4bafc26df99217bf2a1ed1ef45')
 Children:
 1f3450813754300f67abd3d92ae543d1066732de
 Parents:
 a355e47496ff8e6d78a5da6dacb5b38967d28f1a
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/binresol.lib
ra355e47 r555878 2 2 //////////////////////////////////////////////////////////////////////////// 3 3 version="$Id$"; 4 category=" Resolution of singularities";4 category="Commutaive algebra"; 5 5 info=" 6 6 LIBRARY: binresol.lib Combinatorial algorithm of resolution of singularities … … 11 11 @* G. Pfister, pfister@mathematik.unikl.de 12 12 13 MAIN PROCEDURE:14 BINresol(J); computes a Eresolution of singularities of (J) (THE SECOND PART IS NOT IMPLEMENTED YET)15 16 13 PROCEDURES: 17 14 Eresol(J); computes a Eresolution of singularities of (J) in char 0 18 15 BINresol(J); computes a Eresolution of singularities of (J) (THE SECOND PART IS NOT IMPLEMENTED YET) 19 16 determinecenter(L1,L2,c,n,Y,a,mb,flag,control3); computes the next blowingup center 20 21 17 Blowupcenter(L1,id,m,L2,c,n,h); makes the blowingup 22 18 Nonhyp(Coef,expJ,sJ,n,flag,sums); computes the ideal generated by the non hyperbolic generators of expJ 23 24 AUXILIARY PROCEDURES:25 19 inidata(K,k); verifies input data, a binomial ideal K of k generators 26 20 identifyvar(); identifies status of variables 27 28 21 data(K,k,n); transforms data on lists of lenght n 29 22 Edatalist(Coef,Exp,k,n,flag); gives the Eorder of each term in Exp 30 23 EOrdlist(Coef,Exp,k,n,flag); computes the Eorder of an ideal (giving in the language of lists) 31 32 24 maxEord(Coef,Exp,k,n,flag); computes de maximum Eorder of an ideal given by Coef and Exp 33 34 ECoef(Coef,expP,sP,V,auxc,n,flag); Computes a simplified version of the ECoeff ideal. The Eorders are correct, 35 but tranformations of coefficients of the generators and powers of binomials 36 cannot be computed easily in terms of lists. 37 25 ECoef(Coef,expP,sP,V,auxc,n,flag); Computes a simplified version of the ECoeff ideal. 38 26 elimrep(L); removes repeated terms from a list 39 27 Emaxcont(Coef,Exp,k,n,flag); computes a list of hypersurfaces of Emaximal contact 40 28 cleanunit(mon,n,flag); clean the units in a monomial mon 41 42 29 resfunction(t,auxinv,nchart,n); composes the Eresolution function 43 44 30 calculateI(Coef,J,c,n,Y,a,b,D); computes the order of the non monomial part of an ideal J 45 31 Maxord(L,n); computes the maximum exponent of an exceptional monomial ideal 46 32 Gamma(L,c,n); computes the Gamma function for an exceptional monomial ideal given by L 47 48 TRANSLATION OF THE OUTPUT:49 33 convertdata(C,L,n,flag); computes the ideal corresponding to C,L 50 34 tradblwup(blwhist,n,Y,a,num); composes the blowing up at this chart 51 52 35 lcmofall(nchart,mobile); computes the lcm of the denominators of the Eorders for all the charts 53 36 computemcm(Eolist); computes the lcm of the denominators of the Eorders for one chart 54 55 37 constructH(Hhist,n,flag); construct the list of exceptional divisors accumulated at this chart 56 constructblwup(blwhist,n,chy,flag); construct the ideal defining the map K[W] > K[Wi], 57 which gives the composition map of all the blowing up leading to this chart 38 constructblwup(blwhist,n,chy,flag); construct the ideal defining the composition map 58 39 constructlastblwup(blwhist,n,chy,flag); construct the ideal defining the last blowup leading to this chart 59 60 61 40 genoutput(chart,mobile,nchart,nsons,n,q); generates the output for visualization 62 41 salida(idchart,chart,mobile,numson,previousa,n,q); generates the output for one chart 63 64 65 COMPUTATIONS WITH LISTS:66 42 iniD(n); creates a list of lists of zeros of size n 67 43 sumlist(L1,L2); sums two lists component to component … … 71 47 createlist(L1,L2); creates a list of lists of two elements 72 48 list0(n); creates a list of zeros of size n 73 74 49 "; 75 50 … … 377 352 elimrep(L); 378 353 } 379 //////////////////////////////////////////////////////380 381 proc Emaxcont(list Coef,list Exp,int k,int n,list flag)382 "USAGE: Emaxcont(Coef,Exp,k,n,flag);383 Coef,Exp,flag lists, k,n, integers384 Exp is a list of lists of exponents, k=size(Exp)385 COMPUTE: Identify ALL the variables of Emaximal contact386 RETURN: a list with the indexes of the variables of Emaximal contact387 EXAMPLE: example Emaxcont; shows an example388 "389 {390 int i,j,lon;391 number maxEo;392 list L,sums,bx,maxvar;393 394 L=maxEord(Coef,Exp,k,n,flag);395 396 maxEo=L[1];397 sums=L[2];398 399 if (maxEo>0)400 {401 for (i=1;i<=k; i++)402 {403 lon=size(sums[i]);404 if (lon==2)405 {406 if (sums[i][1]==maxEo) // variables of the first term {407 for (j=1;j<=n; j++)408 {409 if(Exp[i][1][j]!=0 and flag[j]==0)410 {411 bx=j; maxvar=maxvar + bx;412 }413 }414 }415 416 if (sums[i][2]==maxEo) // variables of the second term {417 for (j=1;j<=n; j++)418 {419 if(Exp[i][2][j]!=0 and flag[j]==0){bx=j; maxvar=maxvar + bx;}420 }421 }422 }423 else424 {425 if (sums[i][1]==maxEo)426 {427 for (j=1;j<=n; j++)428 {429 if(Exp[i][1][j]!=0 and flag[j]==0)430 {431 bx=j; maxvar=maxvar + bx;432 }433 }434 }435 }436 }437 }438 else {maxvar=list();}439 440 // eliminating repeated terms441 maxvar=elimrep(maxvar);442 443 // It is necessary to check if flag[j]==0 in order to avoid the selection of y variables444 445 return(maxEo,maxvar);446 }447 example448 {"EXAMPLE:"; echo = 2;449 ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;450 list flag=identifyvar();451 ideal J=x(1)^3*x(3)y(2)*y(4)^2,x(5)*y(2)x(7)*y(4)^2,x(6)^2*(1y(4)*y(8)^5),x(7)^4*y(8)^2;452 list L=data(J,4,8);453 list hyp=Emaxcont(L[1],L[2],4,8,flag);454 hyp[1]; // max Eorder=0455 hyp[2]; // There are no hypersurfaces of Emaximal contact456 457 ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;458 list flag=identifyvar();459 ideal J=x(1)^3*x(3)y(2)*y(4)^2*x(3),x(5)*y(2)x(7)*y(4)^2,x(6)^2*(1y(4)*y(8)^5),x(7)^4*y(8)^2;460 list L=data(J,4,8);461 list hyp=Emaxcont(L[1],L[2],4,8,flag);462 hyp[1]; // the Eorder is 1463 hyp[2]; // {x(3)=0},{x(5)=0},{x(7)=0} are hypersurfaces of Emaximal contact464 465 }466 ///////////////////////////////////////////////////////467 354 468 355 proc cleanunit(list mon,int n,list flaglist) … … 707 594 EXAMPLE: example determinecenter; shows an example 708 595 " 709 {int i,j,rstep,l,mm,cont,cont1,cont2,cont3,a4,sI,sP,V,V2,ccase,b,Mindx,tip; 710 number auxc,a1,a2,ex,maxEo,aux; 711 712 list D,H,auxJ; // lists of D_n,D_n1,...,D_1; H_n,H_n1,...,H_1; J_n,J_n1,...,J_1 713 714 list oldOlist,oldC,oldt,oldD,oldH,allH; // information of the previous step 715 716 list Olist,C,t,Dstar,center,expI,expP,newJ,maxset; 717 718 list maxvar,auxlist,aux3,auxD,auxolist,auxdiv,auxaux,L,rs,auxgamma,auxg2,aux1; // auxiliary lists 719 list auxinvlist,newcoef,EL,Ecoaux,Hplus,transH,Hsum,auxset,sumnewH; // auxiliary lists 720 list auxcoefI; 721 722 intvec oldinfobo7,infobo7; 723 int infaux,leh,leh2,leh3; 724 725 tip=listmb[1]; // It is not used in this procedure, it is used to compute the lcm of the denominators 726 oldOlist=listmb[2]; 727 oldC=listmb[3]; 728 oldt=listmb[4]; // t= resolution function 729 oldD=listmb[5]; 730 731 oldH=listmb[6]; 732 allH=listmb[7]; 733 734 oldinfobo7=listmb[8]; // auxiliary intvec, it is used to define BO[7] 735 736 // inicializating lists 737 Olist=list(); 738 C=list(); 739 auxinvlist=list(); 740 741 auxJ[1]=expJ; 742 rstep=n; // we are in dimension rstep 743 auxc=c; 744 cont=1; 745 746 if (Y==0) {D=iniD(n); H=iniD(n); infobo7=1;} // first center, inicializate previous information 747 748 if (Y!=0 and rstep==n) // In dimension n, D'_n is always of this form 749 { auxdiv=list0(n); 750 Dstar[1]=oldD[1]; 751 752 b=size(a); 753 for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[1][i];}}} 754 Dstar[1][Y]=aux; 755 aux=0; 756 757 auxdiv[Y]=oldOlist[1]oldC[1]; 758 D[1]=sumlist(Dstar[1],auxdiv);} // list defining D_n 759 760 // computing strict transforms of the exceptional divisors H 761 762 if (Y!=0){transH=iniD(n); 763 for (i=1;i<=size(oldH);i++){transH[i]=oldH[i]; transH[i][Y]=0;} // Note: size(oldH)<=n 764 allH[Y]=1;} // transform of H=H_nU...UH_1 765 766 // We put here size(oldH) instead of n because maybe we have not 767 // calculated all the dimensions in the previous step 768 769 // STARTING THE LOOP 770 771 while (rstep>=1) 596 { 597 int i,j,rstep,l,mm,cont,cont1,cont2,cont3,a4,sI,sP,V,V2,ccase,b,Mindx,tip; 598 number auxc,a1,a2,ex,maxEo,aux; 599 600 list D,H,auxJ; // lists of D_n,D_n1,...,D_1; H_n,H_n1,...,H_1; J_n,J_n1,...,J_1 601 602 list oldOlist,oldC,oldt,oldD,oldH,allH; // information of the previous step 603 604 list Olist,C,t,Dstar,center,expI,expP,newJ,maxset; 605 606 list maxvar,auxlist,aux3,auxD,auxolist,auxdiv,auxaux,L,rs,auxgamma,auxg2,aux1; // auxiliary lists 607 list auxinvlist,newcoef,EL,Ecoaux,Hplus,transH,Hsum,auxset,sumnewH; // auxiliary lists 608 list auxcoefI; 609 610 intvec oldinfobo7,infobo7; 611 int infaux,leh,leh2,leh3; 612 613 tip=listmb[1]; // It is not used in this procedure, it is used to compute the lcm of the denominators 614 oldOlist=listmb[2]; 615 oldC=listmb[3]; 616 oldt=listmb[4]; // t= resolution function 617 oldD=listmb[5]; 618 619 oldH=listmb[6]; 620 allH=listmb[7]; 621 622 oldinfobo7=listmb[8]; // auxiliary intvec, it is used to define BO[7] 623 624 // inicializating lists 625 Olist=list(); 626 C=list(); 627 auxinvlist=list(); 628 629 auxJ[1]=expJ; 630 rstep=n; // we are in dimension rstep 631 auxc=c; 632 cont=1; 633 634 if (Y==0) {D=iniD(n); H=iniD(n); infobo7=1;} // first center, inicializate previous information 635 636 if (Y!=0 and rstep==n) // In dimension n, D'_n is always of this form 637 { auxdiv=list0(n); 638 Dstar[1]=oldD[1]; 639 640 b=size(a); 641 for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[1][i];}}} 642 Dstar[1][Y]=aux; 643 aux=0; 644 645 auxdiv[Y]=oldOlist[1]oldC[1]; 646 D[1]=sumlist(Dstar[1],auxdiv);} // list defining D_n 647 648 // computing strict transforms of the exceptional divisors H 649 650 if (Y!=0){transH=iniD(n); 651 for (i=1;i<=size(oldH);i++){transH[i]=oldH[i]; transH[i][Y]=0;} // Note: size(oldH)<=n 652 allH[Y]=1;} // transform of H=H_nU...UH_1 653 654 // We put here size(oldH) instead of n because maybe we have not 655 // calculated all the dimensions in the previous step 656 657 // STARTING THE LOOP 658 659 while (rstep>=1) 660 { 661 if (Y!=0 and rstep!=n) // transformation law of D_i for i<n 662 { 663 if (cont!=0) // the resolution function did not drop in higher dimensions 664 { 665 if (oldt[nrstep]==a1/a2 and c==oldC[1] and control3==0) 666 {auxD=list0(n); 667 auxD[Y]=oldOlist[nrstep+1]oldC[nrstep+1]; 668 Dstar[nrstep+1]=oldD[nrstep+1]; 669 670 for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[nrstep+1][i];}}} 671 Dstar[nrstep+1][Y]=aux; 672 aux=0; 673 674 D[nrstep+1]=sumlist(Dstar[nrstep+1],auxD); 675 676 } 677 else 678 {cont=0; 679 for (j=nrstep+1;j<=n; j++){D[j]=list0(n);} 680 } 681 } 682 } 683 684 // Factorizing J=M*I 685 686 cont1=0; 687 for (i=1;i<=n;i++) {if (D[nrstep+1][i]==0) {cont1=cont1+1;}} // if it fails write: listO(n)[i] 688 689 if (cont1==n) {expI=expJ;} // D[nrstep+1]=0 (is a list of zeros) 690 else { 691 for (i=1;i<=size(expJ);i++) 692 {rs[i]=size(Coef[i]); 693 if (rs[i]==2){ aux1=list(); 694 aux1[1]=reslist(expJ[i][1],D[nrstep+1]); 695 aux1[2]=reslist(expJ[i][2],D[nrstep+1]); 696 expI[i]=aux1;} // binomial 697 else {aux1=list(); 698 aux1[1]=reslist(expJ[i][1],D[nrstep+1]); 699 expI[i]=aux1;}} // monomial 700 } 701 702 // NOTE: coeficients of I = coeficients of J, because I and J differ in a monomial 703 704 // Detecting errors, negative exponents in expI 705 706 sI=size(expI); 707 708 for (i=1;i<=sI;i++) 772 709 { 773 if (Y!=0 and rstep!=n) // transformation law of D_i for i<n 710 rs[i]=size(Coef[i]); 711 if (rs[i]==2) 774 712 { 775 if (cont!=0) // the resolution function did not drop in higher dimensions713 for (j=1;j<=2;j++)i 776 714 { 777 if (oldt[nrstep]==a1/a2 and c==oldC[1] and control3==0) 778 {auxD=list0(n); 779 auxD[Y]=oldOlist[nrstep+1]oldC[nrstep+1]; 780 Dstar[nrstep+1]=oldD[nrstep+1]; 781 782 for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[nrstep+1][i];}}} 783 Dstar[nrstep+1][Y]=aux; 784 aux=0; 785 786 D[nrstep+1]=sumlist(Dstar[nrstep+1],auxD); 787 715 for (l=1;l<=n; l++) 716 { 717 if (expI[i][j][l]<0) 718 { 719 ERROR("the BINOMIAL ideal I has negative components"); 720 } 788 721 } 789 else790 {cont=0;791 for (j=nrstep+1;j<=n; j++){D[j]=list0(n);}792 }793 722 } 794 723 } 795 796 // Factorizing J=M*I 797 798 cont1=0; 799 for (i=1;i<=n;i++) {if (D[nrstep+1][i]==0) {cont1=cont1+1;}} // if it fails write: listO(n)[i] 800 801 if (cont1==n) {expI=expJ;} // D[nrstep+1]=0 (is a list of zeros) 802 else { 803 for (i=1;i<=size(expJ);i++) 804 {rs[i]=size(Coef[i]); 805 if (rs[i]==2){ aux1=list(); 806 aux1[1]=reslist(expJ[i][1],D[nrstep+1]); 807 aux1[2]=reslist(expJ[i][2],D[nrstep+1]); 808 expI[i]=aux1;} // binomial 809 else {aux1=list(); 810 aux1[1]=reslist(expJ[i][1],D[nrstep+1]); 811 expI[i]=aux1;}} // monomial 724 else 725 { 726 for (l=1;l<=n; l++) 727 { 728 if (expI[i][1][l]<0) 729 { 730 "the MONOMIAL ideal I has negative components"; 731 "M ideal"; 732 print(D[nrstep+1]); print(expI); 733 "dimension"; print(rstep); 812 734 } 813 814 // NOTE: coeficients of I = coeficients of J, because I and J differ in a monomial 815 816 // Detecting errors, negative exponents in expI 817 818 sI=size(expI); 819 820 for (i=1;i<=sI;i++) 821 {rs[i]=size(Coef[i]); 822 if (rs[i]==2){for (j=1;j<=2;j++){for (l=1;l<=n; l++) 823 {if (expI[i][j][l]<0) {print("ERROR the BINOMIAL ideal I has negative components"); ~;}}}} 824 else {for (l=1;l<=n; l++) 825 {if (expI[i][1][l]<0) {print("ERROR the MONOMIAL ideal I has negative components"); 826 print("M ideal"); print(D[nrstep+1]); print(expI); print("dimension"); print(rstep); ~; }}} 827 } 828 829 // Compute the maximal Eorder of I 830 831 L=maxEord(Coef,expI,sI,n,flag); 832 maxEo=L[1]; // Eorder of I 833 834 // Inicializating information 835 836 auxolist=maxEo; 837 a1=maxEo; 838 a2=auxc; 839 Olist=Olist+auxolist; // list of new maximal Eorders o_n,o_{n1},...o_1 840 aux3=auxc; 841 C=C+aux3; // list of new critical values c=c_{n+1},c_{n},...c_2 842 843 // It is necessary to check if the first coordinate of the invariant has dropped or not 844 // NOTE: By construction, the first coordinate is always 1 !! 845 // It has dropped is equivalent to: CURRENT C<c of the previous step 846 847 // Calculate new H, this is done for every dimension 848 849 if (Y!=0){a4=size(oldt); 850 if (nrstep+1>a4){cont=0; oldt[nrstep+1]=0; } // VERIFICAR!!!! 851 852 if (cont!=0 and oldt[nrstep+1]==a1/a2 and c==oldC[1] and control3==0){H[nrstep+1]=transH[nrstep+1]; 853 854 // we fill now the value for BO[7] 855 if (oldinfobo7[nrstep+1]==1){leh=size(Hhist); 856 infobo7[nrstep+1]=Hhist[leh];} // suitable index !!! 857 else{ infaux=oldinfobo7[nrstep+1]; 858 infobo7[nrstep+1]=infaux;} // the same as the previous step 859 860 } 861 else { 862 if (rstep<n) {sumnewH=list0(n); 863 for (i=1;i<nrstep+1;i++){sumnewH=sumlist(sumnewH,H[i]);} 864 H[nrstep+1]=reslist(allH,sumnewH);} 865 else {H[nrstep+1]=allH;} 866 867 // we fill the value for BO[7] too, we complete it at the end if necessary 868 infobo7[nrstep+1]=1; 869 } 870 } 871 872 // It is necessary to detect the monomial case AFTER inicializate the information 873 // OTHERWISE WE WILL HAVE EMPTY COMPONENTS IN THE RESOLUTION FUNCTION 874 875 // If maxEo=0 but maxo!=0 MONOMIAL CASE (because ESing(J,c) still !=emptyset) 876 // If maxEo=0 and maxo=0 then I=1, (real) monomial case, the same case for us 877 // NOTE THAT IT DOESN'T MATTER IF THERE IS A pTH POWER OF A HYPERBOLIC EQ, THE EORDER IS ZERO ANYWAY 878 879 if (maxEo==0){auxgamma=Gamma(D[nrstep+1],auxc,n); // Gamma gives (maxlist,gamma,center) 880 auxg2=auxgamma[3]; 881 center=center+auxg2; 882 center=elimrep(center); 883 auxinvlist=auxgamma[2]; break;} 884 885 // Calculate P // P=I+M^{o/(co)} with weight o 886 887 if (maxEo>=auxc) {expP=expI; Mindx=0;} // The coefficients also remain constant 888 else {ex=maxEo/(auxcmaxEo); 889 auxlist=list(); 890 Mindx=1; 891 auxlist[1]=multiplylist(D[nrstep+1],ex); // weighted monomial part: D[nrstep+1]^ex; 892 expP=insert(expI,auxlist); // P=I+D[nrstep+1]^ex; 893 auxcoefI=Coef; 894 Coef=insert(Coef,list(1));} // Adding the coefficient for M 895 896 // NOTE: IT IS NECESSARY TO ADD COEFFICIENT 1 TO THE MONOMIAL PART M 897 // Eord(P_i)=Eord(I_i) so to compute the Eorder of P_i we can compute Eord(I_i) 898 899 // Calculate variables of Emaximal contact, ALWAYS WITH RESPECT TO THE IDEAL I !! 900 901 sP=size(expP); // Can be different from size(expI) 902 903 if (Mindx==1){ maxvar=Emaxcont(auxcoefI,expI,sI,n,flag);} 904 else{ maxvar=Emaxcont(Coef,expP,sP,n,flag);} 905 906 auxc=maxvar[1]; // Eorder of P, critical value for the next step, ALSO VALID auxc=maxEo; 907 if (auxc!=maxEo){print("ERROR, the Eorder is not well computed");} 908 909 maxset=maxvar[2]; 910 center=center + maxset; 911 912 // Cleaning the center: eliminating repeated variables 913 914 center=elimrep(center); 915 916 if (rstep==1) {break;} // Induction finished, is not necessary to compute the rest 917 918 // Calculate Hplus=set of non permissible hypersurfaces 919 // RESET Hplus if c has dropped or we have eliminated hyperbolic generators 920 921 // ES NECESARIO PONER CONDICION DE SI INVARIANTE BAJA O NO??? SI BAJA HPLUS NO SE USA... 922 923 if (Y==0 or c<oldC[1] or control3==1) {Hplus=list0(n);} 924 else {Hsum=list0(n); 925 Hplus=allH; 926 for (i=1;i<=nrstep+1;i++){Hsum=sumlist(Hsum,H[i]);} 927 Hplus=reslist(Hplus,Hsum); // CHEQUEAR QUE NO SALEN 1'S 928 } 929 930 // Taking into account variables of maxset outside of Hplus (so inside Hminus) 931 932 if (Y==0 or c<oldC[1] or control3==1){V=maxset[1]; // Hplus=0 so any variable is permissible 933 maxset=delete(maxset,1);} // eliminating this variable V from maxset 934 else{ 935 // If the invariant remains constant V comes from the previous step 936 937 if (cont!=0 and oldt[nrstep+1]==a1/a2 and c==oldC[1]){ 938 if (Mindx==1){ 939 //USING HPLUS 940 // REMIND THAT IN THIS CASE maxset=HYPERSURFACES OF EMAXIMAL CONTACT FOR I, INSTEAD OF P 941 942 cont2=1; 943 cont3=1; 944 auxset=maxset; 945 while (cont2!=0){mm=auxset[1]; 946 if (Hplus[mm]!=0) {auxset=delete(auxset,1); cont3=cont3+1;} 947 // eliminating non permissible variables from maxset 948 else {cont2=0;}} 949 V=maxset[cont3]; // first permissible variable 950 maxset=delete(maxset,cont3); 951 V2=a[nrstep+1]; // can be different from the variable coming from the previous step 952 } 953 954 // 955 else{ V=a[nrstep+1];} 956 } 957 else {V=maxset[1]; // Hplus=0 so any variable is permissible 958 maxset=delete(maxset,1); 959 } 960 961 } 962 963 // Calculate Eco=ECoeff_V(P) where V is a permissible hypersurface which belongs to the center 964 // Eco can have rational exponents 965 966 Ecoaux=ECoef(Coef,expP,sP,V,auxc,n,flag); 967 968 // SPECIAL CASES: BOLD REGULAR CASE 969 // 970 971 if (Ecoaux[3]==1){ // Eco=EMPTY LIST, Eco=0 AS IDEAL 972 aux1[1]=list0(n); 973 newJ[1]=aux1; // monomial with zero entries, newJ=1 as ideal 974 newcoef[1]=list(1); // the new coefficient is only 1 975 auxaux=list(); 976 auxaux[1]=newJ; 977 auxJ=auxJ+auxaux; // auxJ list of ideals J_i 978 auxinvlist=1; 979 break;} 980 981 // 982 // THIS CASE IS NOT GOING TO APPEAR, BUT WE LEAVE IT TO CHECK COMPUTATIONS 983 984 if (Ecoaux[3]==2 or Ecoaux[3]==3){ // Eco=0 LIST, Eco=1 AS IDEAL 985 aux1[1]=list0(n); 986 newJ[1]=aux1; 987 newcoef[1]=list(1); print("Strange case happens"); ~; 988 auxaux=list(); 989 auxaux[1]=newJ; 990 auxJ=auxJ + auxaux; // auxJ list of ideals J_i 991 auxinvlist=1; 992 break;} 993 // 994 // THIS CASE IS NOT GOING TO APPEAR, BUT WE LEAVE IT TO CHECK COMPUTATIONS 995 996 // P=1 THIS CANNOT HAPPEN SINCE P=1 IFF I=1 (or I is equivalent to 1) 997 // and this is the monomial case, already checked 998 999 if (Ecoaux[3]==4){print("ERROR in ECoef"); break;} 1000 // 1001 1002 // If we are here Ecoaux[3]=0, then continue 1003 1004 // Filling the list of "ideals", auxJ=J_n,J_{n1},... 1005 1006 newJ=Ecoaux[1]; 1007 newcoef=Ecoaux[2]; 1008 1009 auxJ=insert(auxJ,newJ,nrstep+1); // newJ is inserted after nrstep+1 position, so in position nrstep+2 1010 1011 // New input for the loop, if we are here newJ is different from 0 1012 1013 expJ=newJ; 1014 Coef=newcoef; 1015 1016 newJ=list(); 1017 expI=list(); 1018 expP=list(); 1019 rstep=rstep1; // print(size(auxJ)); 1020 } 1021 1022 // EXIT LOOP "while" 1023 // we do NOT construct the center as an ideal because WE USE LISTS 1024 1025 t=dividelist(Olist,C); // resolution function t 1026 1027 // Complete the intvec infobo7 if necessary 1028 1029 if (control3==1){infobo7=1;} // We reset the value after clean hyperbolic equations 1030 leh2=size(Olist); 1031 leh3=size(infobo7); 1032 if (leh3<leh2){for (j=leh3+1;j<=leh2; j++){infobo7[j]=1;}} 1033 1034 // Auxiliary list to complete the resolution function in special cases 1035 if (size(auxinvlist)==0) {auxinvlist[1]=0;} 1036 1037 return(center,auxJ,Olist,C,t,D,H,allH,auxinvlist,infobo7); 735 } 736 } 737 } 738 739 // Compute the maximal Eorder of I 740 741 L=maxEord(Coef,expI,sI,n,flag); 742 maxEo=L[1]; // Eorder of I 743 744 // Inicializating information 745 746 auxolist=maxEo; 747 a1=maxEo; 748 a2=auxc; 749 Olist=Olist+auxolist; // list of new maximal Eorders o_n,o_{n1},...o_1 750 aux3=auxc; 751 C=C+aux3; // list of new critical values c=c_{n+1},c_{n},...c_2 752 753 // It is necessary to check if the first coordinate of the invariant has dropped or not 754 // NOTE: By construction, the first coordinate is always 1 !! 755 // It has dropped is equivalent to: CURRENT C<c of the previous step 756 757 // Calculate new H, this is done for every dimension 758 759 if (Y!=0){a4=size(oldt); 760 if (nrstep+1>a4){cont=0; oldt[nrstep+1]=0; } // VERIFICAR!!!! 761 762 if (cont!=0 and oldt[nrstep+1]==a1/a2 and c==oldC[1] and control3==0){H[nrstep+1]=transH[nrstep+1]; 763 764 // we fill now the value for BO[7] 765 if (oldinfobo7[nrstep+1]==1){leh=size(Hhist); 766 infobo7[nrstep+1]=Hhist[leh];} // suitable index !!! 767 else{ infaux=oldinfobo7[nrstep+1]; 768 infobo7[nrstep+1]=infaux;} // the same as the previous step 769 770 } 771 else { 772 if (rstep<n) {sumnewH=list0(n); 773 for (i=1;i<nrstep+1;i++){sumnewH=sumlist(sumnewH,H[i]);} 774 H[nrstep+1]=reslist(allH,sumnewH);} 775 else {H[nrstep+1]=allH;} 776 777 // we fill the value for BO[7] too, we complete it at the end if necessary 778 infobo7[nrstep+1]=1; 779 } 780 } 781 782 // It is necessary to detect the monomial case AFTER inicializate the information 783 // OTHERWISE WE WILL HAVE EMPTY COMPONENTS IN THE RESOLUTION FUNCTION 784 785 // If maxEo=0 but maxo!=0 MONOMIAL CASE (because ESing(J,c) still !=emptyset) 786 // If maxEo=0 and maxo=0 then I=1, (real) monomial case, the same case for us 787 // NOTE THAT IT DOESN'T MATTER IF THERE IS A pTH POWER OF A HYPERBOLIC EQ, THE EORDER IS ZERO ANYWAY 788 789 if (maxEo==0){auxgamma=Gamma(D[nrstep+1],auxc,n); // Gamma gives (maxlist,gamma,center) 790 auxg2=auxgamma[3]; 791 center=center+auxg2; 792 center=elimrep(center); 793 auxinvlist=auxgamma[2]; break;} 794 795 // Calculate P // P=I+M^{o/(co)} with weight o 796 797 if (maxEo>=auxc) {expP=expI; Mindx=0;} // The coefficients also remain constant 798 else {ex=maxEo/(auxcmaxEo); 799 auxlist=list(); 800 Mindx=1; 801 auxlist[1]=multiplylist(D[nrstep+1],ex); // weighted monomial part: D[nrstep+1]^ex; 802 expP=insert(expI,auxlist); // P=I+D[nrstep+1]^ex; 803 auxcoefI=Coef; 804 Coef=insert(Coef,list(1));} // Adding the coefficient for M 805 806 // NOTE: IT IS NECESSARY TO ADD COEFFICIENT 1 TO THE MONOMIAL PART M 807 // Eord(P_i)=Eord(I_i) so to compute the Eorder of P_i we can compute Eord(I_i) 808 809 // Calculate variables of Emaximal contact, ALWAYS WITH RESPECT TO THE IDEAL I !! 810 811 sP=size(expP); // Can be different from size(expI) 812 813 if (Mindx==1){ maxvar=Emaxcont(auxcoefI,expI,sI,n,flag);} 814 else{ maxvar=Emaxcont(Coef,expP,sP,n,flag);} 815 816 auxc=maxvar[1]; // Eorder of P, critical value for the next step, ALSO VALID auxc=maxEo; 817 if (auxc!=maxEo) 818 { 819 ERROR("the Eorder is not well computed"); 820 } 821 822 maxset=maxvar[2]; 823 center=center + maxset; 824 825 // Cleaning the center: eliminating repeated variables 826 827 center=elimrep(center); 828 829 if (rstep==1) {break;} // Induction finished, is not necessary to compute the rest 830 831 // Calculate Hplus=set of non permissible hypersurfaces 832 // RESET Hplus if c has dropped or we have eliminated hyperbolic generators 833 834 // ES NECESARIO PONER CONDICION DE SI INVARIANTE BAJA O NO??? SI BAJA HPLUS NO SE USA... 835 836 if (Y==0 or c<oldC[1] or control3==1) {Hplus=list0(n);} 837 else {Hsum=list0(n); 838 Hplus=allH; 839 for (i=1;i<=nrstep+1;i++){Hsum=sumlist(Hsum,H[i]);} 840 Hplus=reslist(Hplus,Hsum); // CHEQUEAR QUE NO SALEN 1'S 841 } 842 843 // Taking into account variables of maxset outside of Hplus (so inside Hminus) 844 845 if (Y==0 or c<oldC[1] or control3==1){V=maxset[1]; // Hplus=0 so any variable is permissible 846 maxset=delete(maxset,1);} // eliminating this variable V from maxset 847 else{ 848 // If the invariant remains constant V comes from the previous step 849 850 if (cont!=0 and oldt[nrstep+1]==a1/a2 and c==oldC[1]){ 851 if (Mindx==1){ 852 //USING HPLUS 853 // REMIND THAT IN THIS CASE maxset=HYPERSURFACES OF EMAXIMAL CONTACT FOR I, INSTEAD OF P 854 855 cont2=1; 856 cont3=1; 857 auxset=maxset; 858 while (cont2!=0){mm=auxset[1]; 859 if (Hplus[mm]!=0) {auxset=delete(auxset,1); cont3=cont3+1;} 860 // eliminating non permissible variables from maxset 861 else {cont2=0;}} 862 V=maxset[cont3]; // first permissible variable 863 maxset=delete(maxset,cont3); 864 V2=a[nrstep+1]; // can be different from the variable coming from the previous step 865 } 866 867 // 868 else{ V=a[nrstep+1];} 869 } 870 else {V=maxset[1]; // Hplus=0 so any variable is permissible 871 maxset=delete(maxset,1); 872 } 873 874 } 875 876 // Calculate Eco=ECoeff_V(P) where V is a permissible hypersurface which belongs to the center 877 // Eco can have rational exponents 878 879 Ecoaux=ECoef(Coef,expP,sP,V,auxc,n,flag); 880 881 // SPECIAL CASES: BOLD REGULAR CASE 882 // 883 884 if (Ecoaux[3]==1) 885 { // Eco=EMPTY LIST, Eco=0 AS IDEAL 886 aux1[1]=list0(n); 887 newJ[1]=aux1; // monomial with zero entries, newJ=1 as ideal 888 newcoef[1]=list(1); // the new coefficient is only 1 889 auxaux=list(); 890 auxaux[1]=newJ; 891 auxJ=auxJ+auxaux; // auxJ list of ideals J_i 892 auxinvlist=1; 893 break; 894 } 895 // 896 // THIS CASE IS NOT GOING TO APPEAR, BUT WE LEAVE IT TO CHECK COMPUTATIONS 897 if (Ecoaux[3]==2 or Ecoaux[3]==3) 898 { // Eco=0 LIST, Eco=1 AS IDEAL 899 aux1[1]=list0(n); 900 newJ[1]=aux1; 901 newcoef[1]=list(1); print("Strange case happens"); 902 auxaux=list(); 903 auxaux[1]=newJ; 904 auxJ=auxJ + auxaux; // auxJ list of ideals J_i 905 auxinvlist=1; 906 break; 907 } 908 // 909 // THIS CASE IS NOT GOING TO APPEAR, BUT WE LEAVE IT TO CHECK COMPUTATIONS 910 911 // P=1 THIS CANNOT HAPPEN SINCE P=1 IFF I=1 (or I is equivalent to 1) 912 // and this is the monomial case, already checked 913 914 if (Ecoaux[3]==4){ERROR("ERROR in ECoef");} 915 // 916 // If we are here Ecoaux[3]=0, then continue 917 // Filling the list of "ideals", auxJ=J_n,J_{n1},... 918 newJ=Ecoaux[1]; 919 newcoef=Ecoaux[2]; 920 921 auxJ=insert(auxJ,newJ,nrstep+1); // newJ is inserted after nrstep+1 position, so in position nrstep+2 922 923 // New input for the loop, if we are here newJ is different from 0 924 925 expJ=newJ; 926 Coef=newcoef; 927 928 newJ=list(); 929 expI=list(); 930 expP=list(); 931 rstep=rstep1; // print(size(auxJ)); 932 } 933 // EXIT LOOP "while" 934 // we do NOT construct the center as an ideal because WE USE LISTS 935 t=dividelist(Olist,C); // resolution function t 936 // Complete the intvec infobo7 if necessary 937 if (control3==1){infobo7=1;} // We reset the value after clean hyperbolic equations 938 leh2=size(Olist); 939 leh3=size(infobo7); 940 if (leh3<leh2){for (j=leh3+1;j<=leh2; j++){infobo7[j]=1;}} 941 // Auxiliary list to complete the resolution function in special cases 942 if (size(auxinvlist)==0) {auxinvlist[1]=0;} 943 return(center,auxJ,Olist,C,t,D,H,allH,auxinvlist,infobo7); 1038 944 } 1039 945 example … … 1546 1452 oldid=idchart; 1547 1453 1548 if (EordJ<0) {print("ERROR in J in chart"); print(idchart); ~;break;}1454 if (EordJ<0) {print("ERROR in J in chart"); print(idchart); break;} 1549 1455 1550 1456 … … 1552 1458 // CASE J=1, if we reset c, can happen Eord=c=0 1553 1459 1554 // or if there are hyperbolic equations at the beginning!!! A ÑADIR!!!!1460 // or if there are hyperbolic equations at the beginning!!! AÃADIR!!!! 1555 1461 1556 1462 // if (EordJ==0){auxfchart[1]=chart[idchart]; // WE HAVE FINISHED … … 1755 1661 } // closing else for EordI=0 1756 1662 1757 if (EordI<0) {print("ERROR in chart"); print(idchart); ~;break;}1663 if (EordI<0) {print("ERROR in chart"); print(idchart); break;} 1758 1664 1759 1665 // guardar de momento … … 1784 1690 ideal J=x(1)^2x(2)^3; 1785 1691 list L=Eresol(J); 1786 "Please press return after each break point to see the next element of the output list";1787 1692 L[1][1]; // information of the first chart, L[1] list of charts 1788 ~;1789 1693 L[2]; // list of charts with information about sons 1790 ~;1791 1694 L[3]; // invariant, "#" means solved chart 1792 ~;1793 1695 L[4]; // number of charts, 7 in this example 1794 ~;1795 1696 L[5]; // height corresponding to each chart 1796 ~;1797 1697 L[6]; // number of sons 1798 ~;1799 1698 L[7]; // auxiliary invariant 1800 ~;1801 1699 L[8]; // H exceptional divisors and more information 1802 ~;1803 1700 L[9]; // complete resolution function 1804 1701 … … 1809 1706 L[4]; // 16 charts 1810 1707 L[9]; // complete resolution function 1811 ~;1812 1708 1813 1709 "Third example, write L[i] to see the ith component of the list"; … … 1817 1713 L[4]; // 8 charts, rational exponents 1818 1714 L[9]; // complete resolution function 1819 ~;1820 1715 } 1821 1716 … … 2347 2242 setring RR; 2348 2243 BO; 2349 "press return to see next example"; ~;2350 2244 2351 2245 ring r = 0,(x(1..2)),dp; … … 2357 2251 BO; 2358 2252 showBO(BO); 2359 "press return to see next example"; ~;2360 2253 2361 2254 ring r = 0,(x(1..2)),dp; … … 2381 2274 " 2382 2275 { 2383 int idchart,parent; 2384 list auxlist,solvedrings,totalringlist,previousa; 2385 2386 // chart gives: parent,Y,a,expJ,Coef,flag,Hhist,blwhist,path,hipercoef,hiperexp 2387 // mobile gives: tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; NOTE: Eolist=mobile[2]; 2388 2389 idchart=1; 2390 2391 // first loop, construct list previousa 2392 2393 while (idchart<=nchart) 2394 { 2395 if (idchart==1){previousa[1]=chart[2][3];} 2396 2397 else{ 2398 2399 // if there are no sons, the next center is nothing 2400 2401 if (nsons[idchart]==0){previousa[idchart]=0;} 2402 2403 // always fill the parent 2404 2405 parent=chart[idchart][1]; 2406 previousa[parent]=chart[idchart][3]; 2407 } 2408 idchart=idchart+1; 2409 } 2410 2411 // HERE BEGIN THE LOOP 2412 2413 idchart=1; 2414 2415 while (idchart<=nchart) 2416 { 2417 2418 def auxexit=salida(idchart,chart[idchart],mobile[idchart+1],nsons[idchart],previousa[idchart],n,q); 2419 2420 // we add the ring to the list of all rings 2421 2422 auxlist[1]=auxexit; 2423 totalringlist=totalringlist+auxlist; 2424 2425 // if the chart has no sons, add it to the list of final charts 2426 2427 if (nsons[idchart]==0){solvedrings=solvedrings+auxlist;} 2428 2429 auxlist=list(); 2430 kill auxexit; 2431 2432 idchart=idchart+1; 2433 2434 } // EXIT WHILE 2435 2436 return(solvedrings,totalringlist); 2276 int idchart,parent; 2277 list auxlist,solvedrings,totalringlist,previousa; 2278 2279 // chart gives: parent,Y,a,expJ,Coef,flag,Hhist,blwhist,path,hipercoef,hiperexp 2280 // mobile gives: tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; NOTE: Eolist=mobile[2]; 2281 2282 idchart=1; 2283 2284 // first loop, construct list previousa 2285 2286 while (idchart<=nchart) 2287 { 2288 if (idchart==1){previousa[1]=chart[2][3];} 2289 else 2290 { 2291 // if there are no sons, the next center is nothing 2292 if (nsons[idchart]==0){previousa[idchart]=0;} 2293 // always fill the parent 2294 parent=chart[idchart][1]; 2295 previousa[parent]=chart[idchart][3]; 2296 } 2297 idchart=idchart+1; 2298 } 2299 2300 // HERE BEGIN THE LOOP 2301 2302 idchart=1; 2303 2304 while (idchart<=nchart) 2305 { 2306 def auxexit=salida(idchart,chart[idchart],mobile[idchart+1],nsons[idchart],previousa[idchart],n,q); 2307 // we add the ring to the list of all rings 2308 auxlist[1]=auxexit; 2309 totalringlist=totalringlist+auxlist; 2310 // if the chart has no sons, add it to the list of final charts 2311 if (nsons[idchart]==0){solvedrings=solvedrings+auxlist;} 2312 auxlist=list(); 2313 kill auxexit; 2314 idchart=idchart+1; 2315 } // EXIT WHILE 2316 return(solvedrings,totalringlist); 2437 2317 } 2438 2318 example … … 2446 2326 ResTree(B,iden0[1]); // generates the resolution tree 2447 2327 2448 // Use presentTree(B); to see the final charts 2449 // To see the tree type in another shell 2450 // dot Tjpg ResTree.dot o ResTree.jpg 2451 // /usr/bin/X11/xv ResTree.jpg 2452 2328 // Use presentTree(B); to see the final charts 2329 // To see the tree type in another shell 2330 // dot Tjpg ResTree.dot o ResTree.jpg 2331 // /usr/bin/X11/xv ResTree.jpg 2453 2332 } 2454 2333 ///////////////////////////////////////////////////////////////////// … … 2461 2340 " 2462 2341 { 2463 int m,i,aux,mcmchart; 2464 intvec num; 2465 2466 m=size(Eolist); 2467 2468 if (m==1){mcmchart=int(denominator(Eolist[1])); return(mcmchart);} 2469 2470 if (m>1){num=int(denominator(Eolist[1])); 2471 for (i=2;i<=m;i++){aux=int(denominator(Eolist[i])); 2472 num=num,aux; }} 2473 2474 mcmchart=lcm(num); 2475 2476 return(mcmchart); 2342 int m,i,aux,mcmchart; 2343 intvec num; 2344 2345 m=size(Eolist); 2346 2347 if (m==1) 2348 { 2349 mcmchart=int(denominator(Eolist[1])); 2350 return(mcmchart); 2351 } 2352 2353 if (m>1) 2354 { 2355 num=int(denominator(Eolist[1])); 2356 for (i=2;i<=m;i++) 2357 { 2358 aux=int(denominator(Eolist[i])); 2359 num=num,aux; 2360 } 2361 } 2362 mcmchart=lcm(num); 2363 return(mcmchart); 2477 2364 } 2478 2365 example … … 2483 2370 L[8][2][2]; // maximal Eorder at the first chart 2484 2371 computemcm(L[8][2][2]); 2485 2486 2372 } 2487 2373 ///////////////////////////////////////////////////////////////////// … … 2776 2662 EXAMPLE: example createlist; shows an example 2777 2663 " 2778 {int i,k; 2779 list L,aux; 2780 k=size(L1); 2781 if (size(L2)!=k) {return("ERROR en createlist, lists must have the same size");} 2782 L=list0(k); 2783 for (i=1;i<=k;i++) {if (L1[i]!=0) {aux=L1[i],L2[i]; L[i]=aux;} 2784 else {L=delete(L,i);}} 2785 return(L); 2664 { 2665 int i,k; 2666 list L,aux; 2667 k=size(L1); 2668 if (size(L2)!=k) 2669 {ERROR ("createlist: lists must have the same size");} 2670 L=list0(k); 2671 for (i=1;i<=k;i++) 2672 { 2673 if (L1[i]!=0) 2674 { 2675 aux=L1[i],L2[i]; L[i]=aux; 2676 } 2677 else {L=delete(L,i);} 2678 } 2679 return(L); 2786 2680 } 2787 2681 example … … 2797 2691 EXAMPLE: example list0; shows an example 2798 2692 " 2799 {int i; 2800 list L0; 2801 for (i=1;i<=n;i++) {L0[i]=0;} 2802 return(L0); 2693 { 2694 int i; 2695 list L0; 2696 for (i=1;i<=n;i++) {L0[i]=0;} 2697 return(L0); 2803 2698 } 2804 2699 example … … 2807 2702 } 2808 2703 //////////////////////////////////////////////////////////// 2704 2705 proc Emaxcont(list Coef,list Exp,int k,int n,list flag) 2706 "USAGE: Emaxcont(Coef,Exp,k,n,flag); 2707 Coef,Exp,flag lists, k,n, integers 2708 Exp is a list of lists of exponents, k=size(Exp) 2709 COMPUTE: Identify ALL the variables of Emaximal contact 2710 RETURN: a list with the indexes of the variables of Emaximal contact 2711 EXAMPLE: example Emaxcont; shows an example 2712 " 2713 { 2714 int i,j,lon; 2715 number maxEo; 2716 list L,sums,bx,maxvar; 2717 2718 L=maxEord(Coef,Exp,k,n,flag); 2719 2720 maxEo=L[1]; 2721 sums=L[2]; 2722 2723 if (maxEo>0) 2724 { 2725 for (i=1;i<=k; i++) 2726 { 2727 lon=size(sums[i]); 2728 if (lon==2) 2729 { 2730 if (sums[i][1]==maxEo) // variables of the first term 2731 { 2732 for (j=1;j<=n; j++) 2733 { 2734 if(Exp[i][1][j]!=0 and flag[j]==0) 2735 { 2736 bx=j; maxvar=maxvar + bx; 2737 } 2738 } 2739 } 2740 if (sums[i][2]==maxEo) // variables of the second term 2741 { 2742 for (j=1;j<=n; j++) 2743 { 2744 if(Exp[i][2][j]!=0 and flag[j]==0){bx=j; maxvar=maxvar + bx;} 2745 } 2746 } 2747 } 2748 else 2749 { 2750 if (sums[i][1]==maxEo) 2751 { 2752 for (j=1;j<=n; j++) 2753 { 2754 if(Exp[i][1][j]!=0 and flag[j]==0) 2755 { 2756 bx=j; maxvar=maxvar + bx; 2757 } 2758 } 2759 } 2760 } 2761 } 2762 } 2763 else {maxvar=list();} 2764 2765 // eliminating repeated terms 2766 maxvar=elimrep(maxvar); 2767 // It is necessary to check if flag[j]==0 in order to avoid the selection of y variables 2768 return(maxEo,maxvar); 2769 } 2770 example 2771 {"EXAMPLE:"; echo = 2; 2772 ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; 2773 list flag=identifyvar(); 2774 ideal J=x(1)^3*x(3)y(2)*y(4)^2,x(5)*y(2)x(7)*y(4)^2,x(6)^2*(1y(4)*y(8)^5),x(7)^4*y(8)^2; 2775 list L=data(J,4,8); 2776 list hyp=Emaxcont(L[1],L[2],4,8,flag); 2777 hyp[1]; // max Eorder=0 2778 hyp[2]; // There are no hypersurfaces of Emaximal contact 2779 2780 ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp; 2781 list flag=identifyvar(); 2782 ideal J=x(1)^3*x(3)y(2)*y(4)^2*x(3),x(5)*y(2)x(7)*y(4)^2,x(6)^2*(1y(4)*y(8)^5),x(7)^4*y(8)^2; 2783 list L=data(J,4,8); 2784 list hyp=Emaxcont(L[1],L[2],4,8,flag); 2785 hyp[1]; // the Eorder is 1 2786 hyp[2]; // {x(3)=0},{x(5)=0},{x(7)=0} are hypersurfaces of Emaximal contact 2787 } 2788 ///////////////////////////////////////////////////////
Note: See TracChangeset
for help on using the changeset viewer.