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
-
Singular/LIB/binresol.lib (modified) (20 diffs)
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.uni-kl.de 12 12 13 MAIN PROCEDURE:14 BINresol(J); computes a E-resolution of singularities of (J) (THE SECOND PART IS NOT IMPLEMENTED YET)15 16 13 PROCEDURES: 17 14 Eresol(J); computes a E-resolution of singularities of (J) in char 0 18 15 BINresol(J); computes a E-resolution 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 blowing-up center 20 21 17 Blowupcenter(L1,id,m,L2,c,n,h); makes the blowing-up 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 E-order of each term in Exp 30 23 EOrdlist(Coef,Exp,k,n,flag); computes the E-order of an ideal (giving in the language of lists) 31 32 24 maxEord(Coef,Exp,k,n,flag); computes de maximum E-order of an ideal given by Coef and Exp 33 34 ECoef(Coef,expP,sP,V,auxc,n,flag); Computes a simplified version of the E-Coeff ideal. The E-orders 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 E-Coeff 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 E-maximal contact 40 28 cleanunit(mon,n,flag); clean the units in a monomial mon 41 42 29 resfunction(t,auxinv,nchart,n); composes the E-resolution 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 E-orders for all the charts 53 36 computemcm(Eolist); computes the lcm of the denominators of the E-orders 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 E-maximal contact386 RETURN: a list with the indexes of the variables of E-maximal 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*(1-y(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 E-order=0455 hyp[2]; // There are no hypersurfaces of E-maximal 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*(1-y(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 E-order is 1463 hyp[2]; // {x(3)=0},{x(5)=0},{x(7)=0} are hypersurfaces of E-maximal 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_n-1,...,D_1; H_n,H_n-1,...,H_1; J_n,J_n-1,...,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_n-1,...,D_1; H_n,H_n-1,...,H_1; J_n,J_n-1,...,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[n-rstep]==a1/a2 and c==oldC[1] and control3==0) 666 {auxD=list0(n); 667 auxD[Y]=oldOlist[n-rstep+1]-oldC[n-rstep+1]; 668 Dstar[n-rstep+1]=oldD[n-rstep+1]; 669 670 for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[n-rstep+1][i];}}} 671 Dstar[n-rstep+1][Y]=aux; 672 aux=0; 673 674 D[n-rstep+1]=sumlist(Dstar[n-rstep+1],auxD); 675 676 } 677 else 678 {cont=0; 679 for (j=n-rstep+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[n-rstep+1][i]==0) {cont1=cont1+1;}} // if it fails write: listO(n)[i] 688 689 if (cont1==n) {expI=expJ;} // D[n-rstep+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[n-rstep+1]); 695 aux1[2]=reslist(expJ[i][2],D[n-rstep+1]); 696 expI[i]=aux1;} // binomial 697 else {aux1=list(); 698 aux1[1]=reslist(expJ[i][1],D[n-rstep+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[n-rstep]==a1/a2 and c==oldC[1] and control3==0) 778 {auxD=list0(n); 779 auxD[Y]=oldOlist[n-rstep+1]-oldC[n-rstep+1]; 780 Dstar[n-rstep+1]=oldD[n-rstep+1]; 781 782 for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[n-rstep+1][i];}}} 783 Dstar[n-rstep+1][Y]=aux; 784 aux=0; 785 786 D[n-rstep+1]=sumlist(Dstar[n-rstep+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=n-rstep+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[n-rstep+1][i]==0) {cont1=cont1+1;}} // if it fails write: listO(n)[i] 800 801 if (cont1==n) {expI=expJ;} // D[n-rstep+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[n-rstep+1]); 807 aux1[2]=reslist(expJ[i][2],D[n-rstep+1]); 808 expI[i]=aux1;} // binomial 809 else {aux1=list(); 810 aux1[1]=reslist(expJ[i][1],D[n-rstep+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[n-rstep+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[n-rstep+1]); print(expI); print("dimension"); print(rstep); ~; }}} 827 } 828 829 // Compute the maximal E-order of I 830 831 L=maxEord(Coef,expI,sI,n,flag); 832 maxEo=L[1]; // E-order of I 833 834 // Inicializating information 835 836 auxolist=maxEo; 837 a1=maxEo; 838 a2=auxc; 839 Olist=Olist+auxolist; // list of new maximal E-orders o_n,o_{n-1},...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 (n-rstep+1>a4){cont=0; oldt[n-rstep+1]=0; } // VERIFICAR!!!! 851 852 if (cont!=0 and oldt[n-rstep+1]==a1/a2 and c==oldC[1] and control3==0){H[n-rstep+1]=transH[n-rstep+1]; 853 854 // we fill now the value for BO[7] 855 if (oldinfobo7[n-rstep+1]==-1){leh=size(Hhist); 856 infobo7[n-rstep+1]=Hhist[leh];} // suitable index !!! 857 else{ infaux=oldinfobo7[n-rstep+1]; 858 infobo7[n-rstep+1]=infaux;} // the same as the previous step 859 860 } 861 else { 862 if (rstep<n) {sumnewH=list0(n); 863 for (i=1;i<n-rstep+1;i++){sumnewH=sumlist(sumnewH,H[i]);} 864 H[n-rstep+1]=reslist(allH,sumnewH);} 865 else {H[n-rstep+1]=allH;} 866 867 // we fill the value for BO[7] too, we complete it at the end if necessary 868 infobo7[n-rstep+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 E-Sing(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 p-TH POWER OF A HYPERBOLIC EQ, THE E-ORDER IS ZERO ANYWAY 878 879 if (maxEo==0){auxgamma=Gamma(D[n-rstep+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/(c-o)} with weight o 886 887 if (maxEo>=auxc) {expP=expI; Mindx=0;} // The coefficients also remain constant 888 else {ex=maxEo/(auxc-maxEo); 889 auxlist=list(); 890 Mindx=1; 891 auxlist[1]=multiplylist(D[n-rstep+1],ex); // weighted monomial part: D[n-rstep+1]^ex; 892 expP=insert(expI,auxlist); // P=I+D[n-rstep+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 // E-ord(P_i)=E-ord(I_i) so to compute the E-order of P_i we can compute E-ord(I_i) 898 899 // Calculate variables of E-maximal 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]; // E-order of P, critical value for the next step, ALSO VALID auxc=maxEo; 907 if (auxc!=maxEo){print("ERROR, the E-order 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<=n-rstep+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[n-rstep+1]==a1/a2 and c==oldC[1]){ 938 if (Mindx==1){ 939 //----------------------------USING HPLUS---------------------------------------- 940 // REMIND THAT IN THIS CASE maxset=HYPERSURFACES OF E-MAXIMAL 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[n-rstep+1]; // can be different from the variable coming from the previous step 952 } 953 954 //------------------------------------------------------------------------------- 955 else{ V=a[n-rstep+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=E-Coeff_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_{n-1},... 1005 1006 newJ=Ecoaux[1]; 1007 newcoef=Ecoaux[2]; 1008 1009 auxJ=insert(auxJ,newJ,n-rstep+1); // newJ is inserted after n-rstep+1 position, so in position n-rstep+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=rstep-1; // 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 E-order of I 740 741 L=maxEord(Coef,expI,sI,n,flag); 742 maxEo=L[1]; // E-order of I 743 744 // Inicializating information 745 746 auxolist=maxEo; 747 a1=maxEo; 748 a2=auxc; 749 Olist=Olist+auxolist; // list of new maximal E-orders o_n,o_{n-1},...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 (n-rstep+1>a4){cont=0; oldt[n-rstep+1]=0; } // VERIFICAR!!!! 761 762 if (cont!=0 and oldt[n-rstep+1]==a1/a2 and c==oldC[1] and control3==0){H[n-rstep+1]=transH[n-rstep+1]; 763 764 // we fill now the value for BO[7] 765 if (oldinfobo7[n-rstep+1]==-1){leh=size(Hhist); 766 infobo7[n-rstep+1]=Hhist[leh];} // suitable index !!! 767 else{ infaux=oldinfobo7[n-rstep+1]; 768 infobo7[n-rstep+1]=infaux;} // the same as the previous step 769 770 } 771 else { 772 if (rstep<n) {sumnewH=list0(n); 773 for (i=1;i<n-rstep+1;i++){sumnewH=sumlist(sumnewH,H[i]);} 774 H[n-rstep+1]=reslist(allH,sumnewH);} 775 else {H[n-rstep+1]=allH;} 776 777 // we fill the value for BO[7] too, we complete it at the end if necessary 778 infobo7[n-rstep+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 E-Sing(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 p-TH POWER OF A HYPERBOLIC EQ, THE E-ORDER IS ZERO ANYWAY 788 789 if (maxEo==0){auxgamma=Gamma(D[n-rstep+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/(c-o)} with weight o 796 797 if (maxEo>=auxc) {expP=expI; Mindx=0;} // The coefficients also remain constant 798 else {ex=maxEo/(auxc-maxEo); 799 auxlist=list(); 800 Mindx=1; 801 auxlist[1]=multiplylist(D[n-rstep+1],ex); // weighted monomial part: D[n-rstep+1]^ex; 802 expP=insert(expI,auxlist); // P=I+D[n-rstep+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 // E-ord(P_i)=E-ord(I_i) so to compute the E-order of P_i we can compute E-ord(I_i) 808 809 // Calculate variables of E-maximal 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]; // E-order of P, critical value for the next step, ALSO VALID auxc=maxEo; 817 if (auxc!=maxEo) 818 { 819 ERROR("the E-order 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<=n-rstep+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[n-rstep+1]==a1/a2 and c==oldC[1]){ 851 if (Mindx==1){ 852 //----------------------------USING HPLUS---------------------------------------- 853 // REMIND THAT IN THIS CASE maxset=HYPERSURFACES OF E-MAXIMAL 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[n-rstep+1]; // can be different from the variable coming from the previous step 865 } 866 867 //------------------------------------------------------------------------------- 868 else{ V=a[n-rstep+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=E-Coeff_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_{n-1},... 918 newJ=Ecoaux[1]; 919 newcoef=Ecoaux[2]; 920 921 auxJ=insert(auxJ,newJ,n-rstep+1); // newJ is inserted after n-rstep+1 position, so in position n-rstep+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=rstep-1; // 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)^2-x(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 i-th 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 E-order 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 E-maximal contact 2710 RETURN: a list with the indexes of the variables of E-maximal 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*(1-y(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 E-order=0 2778 hyp[2]; // There are no hypersurfaces of E-maximal 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*(1-y(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 E-order is 1 2786 hyp[2]; // {x(3)=0},{x(5)=0},{x(7)=0} are hypersurfaces of E-maximal contact 2787 } 2788 ///////////////////////////////////////////////////////
Note: See TracChangeset
for help on using the changeset viewer.
