Changeset 737960d in git


Ignore:
Timestamp:
Mar 1, 2021, 12:24:11 PM (3 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'd25190065115c859833252500a64cfb7b11e3a50')
Children:
c95beaa99dc4305cb7dbe128571600a256952295
Parents:
31f12ca584f1b01f47ae377d9783a9818653141c
git-author:
Hans Schoenemann <hannes@mathematik.uni-kl.de>2021-03-01 12:24:11+01:00
git-committer:
Hans Schoenemann <hannes@mathematik.uni-kl.de>2021-03-01 12:24:41+01:00
Message:
grobcov.lib N12
File:
1 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/grobcov.lib

    r31f12ca r737960d  
    11//////////////////////////////////////////////////////////////////////////////
    2 version="version grobcov.lib 4.2.0.0  Derc_2020 "; // $Id$
    3            // version N10;  June 2020;
     2version="version grobcov.lib 4.2.0  February_2021 "; // $Id$;
     3           // version N12;  February 2021;
     4
    45info="
    56LIBRARY:  grobcov.lib
     
    120121          hypothesis and thesis in ADGT.
    121122
    122           This version was finished on 26/5/2020,
    123 
     123          The last version N11 (2021) has improved the routines for locus
     124          and allows to determine a parametric locus.
     125
     126          This version was finished on 1/2/2021,
    124127
    125128NOTATIONS: Before calling any routine of the library grobcov,
     
    176179
    177180locus(G); Special routine for determining the geometrical
    178           locus of points verifying given conditions. Given
    179           a parametric ideal J in Q[x1,..,xn][u1,..,um]
    180           with parameters (x1,..,xn) and variables
    181           (u1,..,um), representing the system determining
    182           the n-dimensional locus with tracer point (x1,..,xn)
    183           verifying certain properties,
    184           one can apply locus to the system F, for obtaining the locus.
     181          locus of points verifying given conditions. To use it, the
     182          ring R=(0,a1,..,ap,x1,..xn),(u1,..um,v1..vn),lp;
     183          must be declared, where
     184          (a1,..ap) are parameters (optative),
     185          (x1,..xn) are the variabes of the tracer point,
     186          (u1,..,um) are auxiliary variables,
     187          (v1,..,vn) are the mover variables.
     188          Then the input to locus must be the
     189          parametric ideal F defined in R.
    185190
    186191          locus provides all the components of the locus and
     
    188193           \"Normal\", \"Special\", \"Accumulation\",
    189194          \"Degenerate\".
    190           The mover  are the u-variables.
     195          The mover variables are the last n variables.
    191196          The user can ventually restrict them to a subset of them
    192197          for geometrical reasons but this can change the true
    193198          taxonomy.
     199          locus also allows to determine a parmetric locus
     200          depending on p parameters a1,..ap using then
     201          the option \"numpar\",p.
    194202
    195203locusdg(G); Is a special routine that determines the
     
    203211          Q[x1,..,xn][t1,..,tm] depending on an ideal of
    204212          constraints C in Q[t1,..,tm]. It computes the
    205           locus of the envelop, and determines the
    206           different components as well as its taxonomy:
     213          locus of the envelop, and detemines the
     214          different components as well as their taxonomies:
    207215          \"Normal\", \"Special\", \"Accumulation\",
    208216          \"Degenerate\". (See help for locus).
     
    257265        If H1=1 then H1 is not considered, and analogously for T1.
    258266
    259 ConsLevels(A); Given a list of locally closed sets, constructs the
     267ConsLevels(A); Given a list of locally colsed sets, constructs the
    260268        canonical representation of the levels of A an its complement.
    261269
     
    313321//    Wibmer's Lemma:  19-9-2016
    314322// Final data October 2016
    315 // Updated locus (messages) // Updated locus (messages)
     323// Updated locus (messages)
    316324// Final data Mars 2017
    317325// Release N4: (public)
     
    323331// Release N10: May 2020.
    324332// Updated locus. New determination of the taxonomies
     333// Release N11: February 2021.
     334// Improved the routines for locus. Accept parametric locus as option.
    325335
    326336//*************Auxiliary routines**************
     
    966976        mu=mu*nu;
    967977        rho=qlc[1]*qlm;
     978        //"T_nu="; nu; "mu="; mu; "rho="; rho;
    968979        p=nu*p-rho*F[i];
    969980        r=nu*r;
    970981        for (j=1;j<=size(F);j++){q[j]=nu*q[j];}
    971982        q[i]=q[i]+rho;
     983        //"T_q[i]="; q[i];
    972984        divoc=1;
    973985      }
     
    979991      p=p-lcp*lpp;
    980992    }
     993    //"T_r="; r; "p="; p;
    981994  }
    982995  list res=r,q,mu;
     
    984997}
    985998example
    986 {
    987   "EXAMPLE:"; echo = 2;
     999{ "RXAMPLE:";echo = 2;
     1000  // Division of a polynom by an ideal
     1001
    9881002  if(defined(R)){kill R;}
    9891003  ring R=(0,a,b,c),(x,y),dp;
    9901004  short=0;
     1005
    9911006  // Divisor
    9921007  poly f=(ab-ac)*xy+(ab)*x+(5c);
     1008
    9931009  // Dividends
    994   ideal F=ax+b,cy+a;
     1010  ideal F=ax+b,
     1011       cy+a;
     1012
    9951013  // (Remainder, quotients, factor)
    9961014  def r=pdivi(f,F);
    9971015  r;
     1016
    9981017  // Verifying the division
    9991018  r[3]*f-(r[2][1]*F[1]+r[2][2]*F[2]+r[1]);
     
    11901209};
    11911210example
    1192 {
    1193   "EXAMPLE:"; echo = 2;
    1194   if(defined(R)){kill R;}
    1195   ring R=(0,a,b,c),(x,y),dp;
    1196   short=0;
    1197   // parametric polynom
    1198   poly f=(b^2-1)*x^3*y+(c^2-1)*x*y^2+(c^2*b-b)*x+(a-bc)*y;
    1199   ideal p=c-1;
    1200   ideal q=a-b;
    1201   // normaform of f on V(p) \ V(q)
    1202   pnormalf(f,p,q);
     1211{ "EXAMPLE:"; echo = 2;
     1212
     1213if(defined(R)){kill R;}
     1214ring R=(0,a,b,c),(x,y),dp;
     1215short=0;
     1216
     1217// parametric polynom
     1218poly f=(b^2-1)*x^3*y+(c^2-1)*x*y^2+(c^2*b-b)*x+(a-bc)*y;
     1219// ideals defining V(p)\V(q)
     1220ideal p=c-1;
     1221ideal q=a-b;
     1222
     1223// pnormaform of f on V(p) \ V(q)
     1224pnormalf(f,p,q);
    12031225}
    12041226
     
    14881510        ideals with P included in Q, representing the set
    14891511        V(P) \ V(Q) = V(N) \ V(M)
    1490 KEYWORDS: locally closed set; canonical form
     1512KEYWORDS: locally closed set; canoncial form
    14911513EXAMPLE:  Crep; shows an example"
    14921514{
     
    15181540}
    15191541example
    1520 {
    1521   "EXAMPLE:"; echo = 2;
     1542{ "EXAMPLE:"; echo = 2;
    15221543  short=0;
    15231544  if(defined(R)){kill R;}
     
    15251546  ideal p=a*b;
    15261547  ideal q=a,b-2;
     1548
    15271549  // C-representation of V(p) \ V(q)
    15281550  Crep(p,q);
     
    15881610       Output: [Comp_1, .. , Comp_s ] where
    15891611       Comp_i=[p_i,[p_i1,..,p_is_i]]
    1590 KEYWORDS: locally closed set; canonical form
     1612KEYWORDS: locally closed set; canoncial form
    15911613EXAMPLE:  Prep; shows an example"
    15921614{
     
    16061628}
    16071629example
    1608 {
    1609   "EXAMPLE:"; echo = 2;
     1630{ "EXAMPLE:"; echo = 2;
    16101631  short=0;
    16111632  if(defined(R)){kill R;}
     
    16131634  ideal p=a*b;;
    16141635  ideal q=a,b-1;
     1636
    16151637  // P-representation of V(p) \ V(q)
    16161638  Prep(p,q);
     
    16561678    }
    16571679  }
     1680  //"T_before="; prep;
    16581681  if (size(prep)==0){prep=list(list(ideal(1),list(ideal(1))));}
     1682  //"T_Prep="; prep;
     1683  //def Lout=CompleteA(prep,prep);
     1684  //"T_Lout="; Lout;
    16591685  return(prep);
    16601686}
     
    16811707       P included in Q, representing the
    16821708       set V(P) \ V(Q)
    1683 KEYWORDS: locally closed set; canonical form
     1709KEYWORDS: locally closed set; canoncial form
    16841710EXAMPLE:  PtoCrep; shows an example"
    16851711{
     
    17011727example
    17021728{
    1703  "EXAMPLE:"; echo = 2;
    1704   if(defined(R)){kill R;}
    1705   ring R=0,(a,b,c),lp;
    1706   short=0;
    1707   ideal p=a*(a^2+b^2+c^2-25);
    1708   ideal q=a*(a-3),b-4;
    1709   // C-representaion of V(p) \ V(q)
    1710   def Cr=Crep(p,q);
    1711   Cr;
    1712   // C-representation of V(p) \ V(q)
    1713   def L=Prep(p,q);
    1714   L;
    1715   PtoCrep(L);
     1729echo = 2;
     1730//EXAMPLE:
     1731
     1732if(defined(R)){kill R;}
     1733ring R=0,(a,b,c),lp;
     1734short=0;
     1735
     1736ideal p=a*(a^2+b^2+c^2-25);
     1737ideal q=a*(a-3),b-4;
     1738
     1739// C-representaion of V(p) \ V(q)
     1740def Cr=Crep(p,q);
     1741Cr;
     1742
     1743// P-representation of V(p) \ V(q)
     1744def L=Prep(p,q);
     1745L;
     1746
     1747PtoCrep(L);
    17161748}
    17171749
     
    17341766  {
    17351767    option(returnSB);
     1768    //"T_Lp[i]="; Lp[i];
    17361769    N=Lp[i][1];
    17371770    Lb=Lp[i][2];
     1771    //"T_Lb="; Lb;
    17381772    ida=intersect(ida,N);
    17391773    for(j=1;j<=size(Lb);j++)
     
    19471981    def BH=imap(RR,B2);
    19481982    def LH=imap(RR,LL);
     1983    //"T_BH="; BH;
     1984    //"T_LH="; LH;
    19491985    for (i=1;i<=size(BH);i++)
    19501986    {
     
    19581994    def RPH=DH+PH;
    19591995    def GSH=KSW(BH,LH);
     1996    //"T_GSH="; GSH;
    19601997    //setglobalrings();
    19611998    // DEHOMOGENIZING THE RESULT
     
    20112048example
    20122049{
    2013   "EXAMPLE:"; echo = 2;
    2014   // Casas conjecture for degree 4:
    2015   if(defined(R)){kill R;}
    2016   ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
    2017   short=0;
    2018   ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0),
    2019           x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1),
    2020           x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0),
    2021           x2^2+(2*a3)*x2+(a2),
    2022           x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0),
    2023           x3+(a3);
    2024   cgsdr(F);
     2050echo = 2;
     2051// EXAMPLE:
     2052// Casas conjecture for degree 4:
     2053
     2054// Casas-Alvero conjecture states that on a field of characteristic 0,
     2055// if a polynomial of degree n in x has a common root whith each of its
     2056// n-1 derivatives (not assumed to be the same), then it is of the form
     2057// P(x) = k(x + a)^n, i.e. the common roots must all be the same.
     2058
     2059if(defined(R)){kill R;}
     2060ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
     2061short=0;
     2062
     2063ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0),
     2064         x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1),
     2065         x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0),
     2066         x2^2+(2*a3)*x2+(a2),
     2067         x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0),
     2068         x3+(a3);
     2069
     2070cgsdr(F);
    20252071}
    20262072
     
    20912137    P[i]=L[Q[i][1]][Q[i][2]];
    20922138  }
     2139  //"T_P="; P;
    20932140  // P is the list of the maximal components of the union
    20942141  //   with the corresponding initial holes.
     
    21282175  int i; int j; int k; list H; list C; list T;
    21292176  list L0; list P0; list P; list Q0; list Q;
     2177  //"T_L="; L;
    21302178  for (i=1;i<=size(L);i++)
    21312179  {
     
    21362184    }
    21372185  }
     2186  //"T_P0="; P0;
    21382187  Q0=selectminideals(P0);
     2188  //"T_Q0="; Q0;
    21392189  for (i=1;i<=size(Q0);i++)
    21402190  {
     
    21422192    P[i]=L[Q0[i][1]];// [Q[i][2]];
    21432193  }
     2194  //"T_P="; P;
    21442195  // P is the list of the maximal components of the union
    21452196  //   with the corresponding initial holes.
     
    22492300static proc addpartfine(list H, list C0)
    22502301{
     2302  //"T_H="; H;
    22512303  int i; int j; int k; int te; intvec notQ; int l; list sel;
    22522304  intvec jtesC;
     
    23352387static proc needcombine(list LCU,ideal N)
    23362388{
     2389  //"Deciding if combine is needed";;
    23372390  ideal BB;
    23382391  int tes=0; int m=1; int j; int k; poly sp;
     
    25092562      LCU[k][1]=B;
    25102563    }
     2564    //"Deciding if combine is needed";
    25112565    crep=PtoCrep(prep);
    25122566    if(size(LCU)>1){tes=1;}
     
    27892843example
    27902844{
    2791 "EXAMPLE:"; echo = 2;
     2845echo = 2;
     2846// EXAMPLE 1:
    27922847// Casas conjecture for degree 4:
    2793   if(defined(R)){kill R;}
    2794   ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
    2795   short=0;
    2796   ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0),
    2797             x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1),
    2798             x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0),
    2799             x2^2+(2*a3)*x2+(a2),
    2800             x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0),
    2801             x3+(a3);
    2802   grobcov(F);
    2803 
    2804 // EXAMPLE:
     2848
     2849// Casas-Alvero conjecture states that on a field of characteristic 0,
     2850// if a polynomial of degree n in x has a common root whith each of its
     2851// n-1 derivatives (not assumed to be the same), then it is of the form
     2852// P(x) = k(x + a)^n, i.e. the common roots must all be the same.
     2853
     2854if(defined(R)){kill R;}
     2855ring R=(0,a0,a1,a2,a3,a4),(x1,x2,x3),dp;
     2856short=0;
     2857
     2858ideal F=x1^4+(4*a3)*x1^3+(6*a2)*x1^2+(4*a1)*x1+(a0),
     2859         x1^3+(3*a3)*x1^2+(3*a2)*x1+(a1),
     2860         x2^4+(4*a3)*x2^3+(6*a2)*x2^2+(4*a1)*x2+(a0),
     2861         x2^2+(2*a3)*x2+(a2),
     2862         x3^4+(4*a3)*x3^3+(6*a2)*x3^2+(4*a1)*x3+(a0),
     2863         x3+(a3);
     2864
     2865grobcov(F);
     2866
     2867// EXAMPLE 2
    28052868// M. Rychlik robot;
    28062869// Complexity and Applications of Parametric Algorithms of
    2807 //    Computational Algebraic Geometry.;
     2870// Computational Algebraic Geometry.;
    28082871// In: Dynamics of Algorithms, R. de la Llave, L. Petzold and J. Lorenz eds.;
    28092872// IMA Volumes in Mathematics and its Applications,
    2810 //    Springer-Verlag 118: 1-29 (2000).;
    2811 //    (18. Mathematical robotics: Problem 4, two-arm robot)."
    2812   if (defined(R)){kill R;}
    2813   ring R=(0,a,b,l2,l3),(c3,s3,c1,s1), dp;
    2814   short=0;
    2815   ideal S12=a-l3*c3-l2*c1,b-l3*s3-l2*s1,c1^2+s1^2-1,c3^2+s3^2-1;
    2816   S12;
    2817   grobcov(S12);
     2873// Springer-Verlag 118: 1-29 (2000).;
     2874// (18. Mathematical robotics: Problem 4, two-arm robot).
     2875
     2876if (defined(R)){kill R;}
     2877ring R=(0,a,b,l2,l3),(c3,s3,c1,s1), dp;
     2878short=0;
     2879
     2880ideal S12=a-l3*c3-l2*c1,b-l3*s3-l2*s1,c1^2+s1^2-1,c3^2+s3^2-1;
     2881S12;
     2882
     2883grobcov(S12);
    28182884}
    28192885
     
    29162982example
    29172983{
    2918   echo = 2; "EXAMPLE:";
    2919   if(defined(R)){kill R;}
    2920   ring R=(0,a1,a2),(x),lp;
    2921   short=0;
    2922   poly f=(a1^2-4*a1+a2^2-a2)*x+(a1^4-16*a1+a2^3-4*a2);
    2923   ideal p=a1*a2;
    2924   ideal q=a2^2-a2,a1*a2,a1^2-4*a1;
    2925   extendpoly(f,p,q);" ";
    2926 // EXAMPLE;
     2984echo = 2;
     2985// EXAMPLE 1
     2986
     2987if(defined(R)){kill R;}
     2988ring R=(0,a1,a2),(x),lp;
     2989short=0;
     2990
     2991poly f=(a1^2-4*a1+a2^2-a2)*x+(a1^4-16*a1+a2^3-4*a2);
     2992ideal p=a1*a2;
     2993ideal q=a2^2-a2,a1*a2,a1^2-4*a1;
     2994
     2995extendpoly(f,p,q);
     2996
     2997// EXAMPLE 2
     2998
    29272999if (defined(R)){kill R;}
    29283000ring R=(0,a0,b0,c0,a1,b1,c1,a2,b2,c2),(x), dp;
    29293001short=0;
     3002
    29303003poly f=(b1*a2*c2-c1*a2*b2)*x+(-a1*c2^2+b1*b2*c2+c1*a2*c2-c1*b2^2);
    29313004ideal p=
     
    30903163example
    30913164{
    3092   "EXAMPLE:"; echo = 2;
     3165echo = 2;
     3166// EXAMPLE
     3167
    30933168if(defined(R)){kill R;}
    30943169ring R=(0,a0,b0,c0,a1,b1,c1),(x), dp;
    30953170short=0;
     3171
    30963172ideal S=a0*x^2+b0*x+c0,
    30973173          a1*x^2+b1*x+c1;
     3174
    30983175def GCS=grobcov(S,"rep",2);
    3099   // grobcov(S) with both P and C representations
     3176// grobcov(S) with both P and C representations
    31003177GCS;
    31013178
    31023179def FGC=extendGC(GCS,"rep",1);
    3103   // Full representation
     3180// Full representation
    31043181FGC;
    31053182}
     
    32433320  }
    32443321  combpolys=reform(combpolys,numdens);
     3322  //"T_combpolys="; combpolys;
    32453323  //setring(RR);
    32463324  //def combpolysT=imap(P,combpolys);
     
    33763454//    where the w_l=(n_l1,..,n_ls) are intvec of length size(L), where
    33773455//    n_lt fixes which element of (v_t1,..,v_tk_t) is to be
    3378 //    chosen to form the tth (Q,P) for the lth element of the sheaf
     3456//    choosen to form the tth (Q,P) for the lth element of the sheaf
    33793457//    representing the I-regular function.
    33803458// The selection is done to obtian the minimal number of elements
     
    39564034  {
    39574035    CG=KSWtocgsdr(CG);
     4036    //"T_CG="; CG;
    39584037    if( size(CG)>0)
    39594038    {
     
    43494428  {
    43504429    t=T[i];
     4430    //"T_t"; t;
    43514431    P=0; Q=1;
    43524432    for(j=1;j<=n;j++)
     
    43624442    }
    43634443    Cb=Crep0(P,Q);
     4444    //"T_Cb="; Cb;
    43644445    if(size(Cb)!=0)
    43654446    {
     
    44504531  {
    44514532    K=FirstLevel(B);
     4533    //"T_K="; K;
    44524534    L[size(L)+1]=K[1];
    44534535    B=K[2];
     
    44604542example
    44614543{
    4462 "EXAMPLE:"; echo = 2;
     4544echo = 2;
     4545// EXAMPLE:
     4546
    44634547if(defined(R)){kill R;}
    44644548ring R=0,(x,y,z),lp;
    44654549short=0;
     4550
    44664551ideal P1=(x^2+y^2+z^2-1);
    44674552ideal Q1=z,x^2+y^2-1;
     
    44794564def LL=ConsLevels(L);
    44804565LL;
     4566
    44814567def LLL=Levels(LL);
    44824568LLL;
     
    45174603example
    45184604{
    4519 "EXAMPLE:"; echo = 2;
     4605echo = 2;
     4606// EXAMPLE:
     4607
    45204608if(defined(R)){kill R;}
    45214609ring R=0,(x,y,z),lp;
    45224610short=0;
     4611
    45234612ideal P1=(x^2+y^2+z^2-1);
    45244613ideal Q1=z,x^2+y^2-1;
     
    45364625def LL=ConsLevels(L);
    45374626LL;
     4627
    45384628def LLL=Levels(LL);
    45394629LLL;
     
    45704660  if (size(B) mod 2==1){B[size(B)+1]=ideal(1);}
    45714661  if (size(A) mod 2==1){A[size(A)+1]=ideal(1);}
     4662  //"T_A=";A;
     4663 // "T_B="; B;
    45724664  n=size(A) div 2;
    45734665  m=(size(B) div 2)-1;
     
    45884680      //string("T_j=",j);
    45894681      ABdw=intersectpar(list(A[2*i],B[2*j+1]));
     4682      //"T_ABdw="; ABdw;
    45904683      ABup=A[2*i-1];
     4684      //"T_ABup1="; ABup;
    45914685      if(j>0)
    45924686      {
     
    45964690        }
    45974691      }
     4692      //"T_ABup2="; ABup;
    45984693      ABupC=Crep(ABup,ideal(1));
     4694      //"T_ABupC="; ABupC;
    45994695      ABup=ABupC[1];
     4696       //"T_ABup="; ABup;
    46004697      if(ABup==1){t=0;}
    46014698     //if(equalideals(ABup,ideal(1))){t=0;}
     
    46034700      {
    46044701        M=Crep(ABup,ABdw);
     4702        //"T_M="; M;
    46054703        //if(not(equalideals(M[1],ideal(1)))) {L[size(L)+1]=M;}
    46064704        if(not(size(M)==0)) {L[size(L)+1]=M;}
    46074705      }
     4706      //"L="; L;
    46084707      j++;
    46094708    }
     
    46144713example
    46154714{
    4616 "EXAMPLE:"; echo = 2;
     4715echo = 2;
     4716// EXAMPLE:
     4717
    46174718if(defined(R)){kill R;}
    46184719ring R=(0,x,y,z,t),(x1,y1),lp;
     
    46354736def LL=DifConsLCSets(L1,L2);
    46364737LL;
     4738
    46374739def LLL=ConsLevels(LL);
    46384740LLL;
     4741
    46394742def LLLL=Levels(LLL);
    46404743LLLL;
     
    46444747
    46454748//******************** Begin locus and envelop ******************************
     4749
     4750// Routines for defining different rings acting in the basic ring RR=Q[a,x][u,v], in lp order, where
     4751// a= parameters of the locus problem
     4752// x= tracer variables
     4753// u= auxiiary variables
     4754// v= mover variables
     4755
     4756// Transforms the ringlist of Q[x_1,..,x_j] into the ringlist of Q[x_1,..,x_{n-1},x_{m+1},..,x_j]
     4757//  I.e., deletes the varibles x_n  to x_m
     4758// To be used with the  same order for all variables
     4759static proc  Ldelnm(list LQx,int  n,int  m)
     4760{
     4761  int i;
     4762  int npara= m- n+1;
     4763  def RR=basering;
     4764  def LR=LQx;
     4765  int nt=size(LR[2]);
     4766  def L1=LR[2];
     4767  for(i=n;i<= m;i++) {L1=delete(L1,n);}
     4768  LR[2]=L1;
     4769  intvec v;
     4770  for(i=1;i<=nt-npara;i++){v[i]=1;}
     4771  LR[3][1][2]=v;
     4772  return(LR);
     4773}
     4774
     4775// Transforms the ringlist of Q[a][x] into the ringlist of Q[a,x]
     4776// To be used with the same lp order
     4777proc La_xToLax(list La_x)
     4778{
     4779  if(typeof(La_x[1])==typeof(0)){return(La_x);}
     4780  list Lax=La_x[1];
     4781  if(Lax[1]=0){return(La_x);}
     4782  list Va=Lax[2];
     4783  int na=size(Lax[2]);
     4784  //"na=";na;
     4785  list Vx=La_x[2];
     4786  list Vax=Va+Vx;
     4787  int nx=size(Vx);
     4788  //"nx="; nx;
     4789  intvec vv;
     4790  int i;
     4791  for(i=1;i<=na+nx;i++){vv[i]=1;}
     4792  Lax[2]=Vax;
     4793  Lax[3][1][2]=vv;
     4794  return(Lax);
     4795}
     4796
     4797// Transforms the ringlist of Q[a,x] into the ringlist of Q[a][x]
     4798// To be used with the same lp order
     4799proc LaxToLa_x(list Lax,int nx)
     4800{
     4801  //"T_Lax=",Lax;
     4802  int nax=size(Lax[2]);
     4803  int na=nax-nx;
     4804  if(na==0){return(Lax);}
     4805  else
     4806  {
     4807    //string("T_ na=",na,", nx=",nx);
     4808    list La_x;
     4809    list Vax=Lax[2];
     4810    list Va;
     4811    list Vx;
     4812    int i;
     4813    for(i=1;i<=na;i++){Va[size(Va)+1]=Vax[i];}
     4814    intvec vva;
     4815    for(i=1;i<=na;i++){vva[i]=1;}
     4816    intvec vvx;
     4817    for(i=1;i<=nx;i++){Vx[size(Vx)+1]=Vax[na+i];}
     4818    for(i=1;i<=nx;i++){vvx[i]=1;}
     4819    La_x[1]=Lax;
     4820    La_x[1][2]=Va;
     4821    La_x[1][3][1][2]=vva;
     4822    La_x[2]=Vx;
     4823    list lax3;
     4824    lax3=Lax[3];
     4825    lax3[1][2]=vvx;
     4826    La_x[3]=lax3;
     4827    La_x[4]=Lax[4];
     4828    return(La_x);
     4829  }
     4830}
     4831
     4832//  // Transforms the ringlist of Q[a,x] into the ringlist of Q[a][x]
     4833//  // To be used with the same lp order
     4834//  proc LaxToLa_x(list Lax,int nx)
     4835//  {
     4836//    //"T_Lax=",Lax;
     4837//    int nax=size(Lax[2]);
     4838//    int na=nax-nx;
     4839//    if(na==0){return(Lax);}
     4840//    else
     4841//    {
     4842//      //string("T_ na=",na,", nx=",nx);
     4843//      list La_x;
     4844//      list Vax=Lax[2];
     4845//      list Va;
     4846//      list Vx;
     4847//      int i;
     4848//      for(i=1;i<=na;i++){Va[size(Va)+1]=Vax[i];}
     4849//      intvec vva;
     4850//      for(i=1;i<=na;i++){vva[i]=1;}
     4851//      intvec vvx;
     4852//      for(i=1;i<=nx;i++){Vx[size(Vx)+1]=Vax[na+i];}
     4853//      for(i=1;i<=nx;i++){vvx[i]=1;}
     4854//      La_x[1]=Lax;
     4855//      La_x[1][2]=Va;
     4856//      La_x[1][3][1][2]=vva;
     4857//      La_x[2]=Vx;
     4858//      list lax3;
     4859//      lax3=Lax[3];
     4860//      lax3[1][2]=vvx;
     4861//      La_x[3]=lax3;
     4862//      La_x[4]=Lax[4];
     4863//      return(La_x);
     4864//    }
     4865//  }
     4866
     4867//  proc Lax_uvToLa_v(list Lax_uv,int na, int nv)
     4868//  {
     4869//    int i;
     4870//    def Lax=Lax_uv[1];
     4871//    int nax=size(Lax[2]);
     4872//    int nuv=size(Lax_uv[2]);
     4873//    def La=Lax;
     4874//    list La2;
     4875//    if(na===0){
     4876//      list Luv=Lax_uv;
     4877//      Luv[1]=0;
     4878//      list Lv=
     4879//      ;
     4880//
     4881//    }
     4882//    for(i=1;i<=na;i++){}
     4883//  }
     4884
     4885// Transforms the set of ringlists of Q[a] and Q[x] into the ringlist of Q[a][x]
     4886// To be used with the same lp order
     4887proc LaLxToLa_x(list La,list Lx)
     4888{
     4889  if(size(La)==0){return(Lx);}
     4890  def L1=La;
     4891  def L2=Lx;
     4892  list L;
     4893  L[1]=L1;
     4894  L[2]=L2[2];
     4895  L[3]=L2[3];
     4896  L[4]=L2[4];
     4897  return(L);
     4898}
     4899
     4900// Transforms the ringlist of Q[a,x] into the ring of Q[a]
     4901//  proc LaxToLa(list Lax, int na)
     4902//  {
     4903//    int ntot=size(Lax[2]);
     4904//    list La=Lax;
     4905//    int i;
     4906//    list V;
     4907//    for(i=1;i<=ntot-na;i++){V[i]=Lax[2][na+i];}
     4908//    La[2]=V;
     4909//    intvec vv;
     4910//    for(i=1;i<=ntot-na;i++){vv[i]=1;}
     4911//    La[3][1][2]=vv;
     4912//    return(La);
     4913//  }
     4914
     4915
     4916// Transforms the set of ringlists of Q[a] and Q[x] into the ringlist of Q[a,x]
     4917// To be used with the same lp order
     4918proc LaLxToLax(list La,list Lx)
     4919{
     4920  list L=LaLxToLa_x(La,Lx);
     4921  list Lax=La_xToLax(L);
     4922  return(Lax);
     4923}
     4924
     4925// Transforms the ringlist of Q[a,x] into the ringlist of Q[a]
     4926// To be used with the same lp order
     4927proc LaxToLa(list Lax,int na)
     4928{
     4929  list La;
     4930  if(na==0){return(La);}
     4931  else
     4932  {
     4933    La=Lax;
     4934    list La2;
     4935    for(i=1;i<=na;i++){La2[i]=Lax[2][i];}
     4936    La[2]=La2;
     4937    intvec va;
     4938    for(i=1;i<=na;i++){va[i]=1;}
     4939    La[3][1][2]=va;
     4940    return(La);
     4941  }
     4942}
     4943
     4944// Transforms the ringlist of Q[a,x][u,v] into the ringlist of Q[a][v]
     4945// To be used with the same lp order
     4946proc Lax_uvToLa_v(list Lax_uv,int na, int nv)
     4947{
     4948  //string("T_na=",na,"; nv=",nv);
     4949  int i;
     4950  list Lax=Lax_uv[1];
     4951  int nax=size(Lax[2]);
     4952  int nuv=size(Lax_uv[2]);
     4953  list Lv=Lax_uv;
     4954  int nx=nax-na;
     4955  int nu=nuv-nv;
     4956  Lv[1]=0;
     4957  list Lv2;
     4958  intvec vv;
     4959  for(i=1;i<=nv;i++){Lv2[i]=Lv[2][nu+i];}
     4960  for(i=1;i<=nv;i++){vv[nu+i]=1;}
     4961  Lv[2]=Lv2;
     4962  Lv[3][1][2]=vv;
     4963  if(na==0){return(Lv);}
     4964  else
     4965  {
     4966    list La=Lax;
     4967    list La2;
     4968    intvec va;
     4969    for(i=1;i<=na;i++){La2[i]=Lax[2][i];}
     4970    for(i=1;i<=na; i++){va[i]=1;}
     4971    La[2]=La2;
     4972    La[3][1][2]=va;
     4973    //"T_La="; La;
     4974    //"T_Lv="; Lv;
     4975    list La_v=LaLxToLa_x(La,Lv);
     4976    return(La_v);
     4977  }
     4978}
     4979
     4980// Transforms the ringlist of Q[a,x] [u,v]into the ringlist of Q[x,u,a,v]
     4981// To be used with the same lp order
     4982proc Lax_uvToLxuav(list Lax_uv, int na, int nv)
     4983{
     4984  //string("T_na=",na,"; nv=",nv);
     4985  int i;
     4986  //"T_Lax_uv="; Lax_uv;
     4987  int nax=size(Lax_uv[1][2]);
     4988  int nuv=size(Lax_uv[2]);
     4989  int nx=nax-na;
     4990  int nu=nuv-nv;
     4991  //string("T_nax=",nax,"; nuv=",nuv,"; nx=",nx,"; nu=",nu);
     4992  list Lxuav=Lax_uv[1];
     4993  list L2;
     4994  for(i=1;i<=nx;i++){L2[i]=Lax_uv[1][2][na+i];}
     4995  for(i=1;i<=nu;i++){L2[nx+i]=Lax_uv[2][i];}
     4996  for(i=1;i<=na;i++){L2[nx+nu+i]=Lax_uv[1][2][i];}
     4997  for(i=1;i<=nv;i++){L2[nx+nu+na+i]=Lax_uv[2][nu+i];}
     4998  Lxuav[2]=L2;
     4999  intvec vv;
     5000  for(i=1;i<=nax+nuv;i++){vv[i]=1;}
     5001  Lxuav[3][1][2]=vv;
     5002  return(Lxuav);
     5003}
     5004
     5005
    46465006
    46475007// indepparameters
     
    47195079  attrib(NP,"IsSB",1);
    47205080  int d=dim(std(NP));
     5081  //"T_d="; d;
    47215082  if(d==0){te=0;}
    47225083  setring(RR);
     
    47245085}
    47255086
    4726 //  DimPar
    4727 //  Auxiliary routine to NorSing determining the dimension of a parametric ideal
    4728 //  Does not use @P and define it directly because is assumes that
    4729 //                              variables and parameters have been inverted
    4730  static proc DimPar(ideal E)
     5087//  DimPar(E,nax,nx):
     5088//  Auxilliary routine of locus2 determining the dimension of a component of the locus in
     5089//  the ring Q[a][x]
     5090 static proc DimPar(ideal E,nax,nx)
    47315091 {
     5092   //" ";"T_E in DimPar="; E;
    47325093   def RRH=basering;
    47335094   def RHx=ringlist(RRH);
    47345095   def P=ring(RHx[1]);
    4735    setring(P);
     5096   list Lax=ringlist(P);
     5097   //"Lax="; Lax;
     5098   //int nax=size(Lax[2]);
     5099   int na=nax-nx;
     5100   //string("T_na=",na,"; nx=",nx);
     5101   list La_x=LaxToLa_x(Lax,nx);
     5102   //"T_La_x="; La_x;
     5103   def Qa_x=ring(La_x);
     5104   setring(Qa_x);
     5105   //setring(P);
    47365106   def E2=std(imap(RRH,E));
     5107   //"T_E2 in DimPar="; E2;
    47375108   attrib(E2,"IsSB",1);
    47385109   def d=dim(E2);
     5110   //string("T_d in DimPar=", d);" ";
    47395111   setring RRH;
    47405112   return(d);
    47415113 }
    47425114
    4743 //  DimComp
    4744 //  Auxiliary routine to locus2 determining the dimension of a parametric ideal
    4745 //  it is identical to DimPar but adds information about the character of the component
    4746 static proc dimComp(ideal PA)
    4747 {
    4748   def RR=basering;
    4749   list Rx=ringlist(RR);
    4750   int n=size(Rx[1][2]);
    4751   def P=ring(Rx[1]);
    4752   setring(P);
    4753   list Lout;
    4754   def CP=imap(RR,PA);
    4755   attrib(CP,"IsSB",1);
    4756   int d=dim(std(CP));
    4757   if(d==n-1){Lout[1]=d;Lout[2]="Degenerate"; }
    4758   else {Lout[1]=d; Lout[2]="Accumulation";}
    4759   setring RR;
    4760   return(Lout);
     5115//     DimComp
     5116//     Auxilliary routine of locus2 determining the dimension of a parametric ideal
     5117//     it is identical to DimPar but adds infromation about the character of the component
     5118 static proc DimComp(ideal PA, int nax,int nx)
     5119 {
     5120//     def RR=basering;
     5121//     list Rx=ringlist(RR);
     5122//     int nax=size(Rx[1][2]);
     5123//     int na=nax-nx;
     5124//     def P=ring(Rx[1]);
     5125//     setring(P);
     5126//     list Lout;
     5127//     def CP=imap(RR,PA);
     5128//     attrib(CP,"IsSB",1);
     5129//     int d=dim(std(CP));
     5130
     5131   list Lout;
     5132   int d=DimPar(PA,nax,nx);
     5133   if(d==nax-1){Lout[1]=d;Lout[2]="Degenerate"; }
     5134   else {Lout[1]=d; Lout[2]="Accumulation";}
     5135   //"T_Lout="; Lout;
     5136   setring RR;
     5137   return(Lout);
    47615138}
    47625139
     
    47825159  return(L3);
    47835160}
    4784 
    4785 // NorSing
    4786 // Input:
    4787 //   ideal B: the basis of a component of the grobcov
    4788 //   ideal E: the top of the component (assumed to be of dimension > 0 (single equation)
    4789 //   ideal N: the holes of the component
    4790 // Output:
    4791 //   int d: the dimension of B on the variables (anti-image).
    4792 //     if d>0    then the component is 'Normal'
    4793 //     if d==0 then the component is 'Singular'
    4794 static proc NorSing(ideal B, ideal E, ideal N, list #)
    4795   {
    4796     int i; int j; int Fenv=0; int env; int dd;
    4797     list DD=#;
    4798     def RR=basering;
    4799     ideal vmov;
    4800     int version=0;
    4801     int nv=nvars(RR);
    4802     int moverdim=nv; //npars(RR);
    4803     if(nv<4){version=1;}
    4804     int d;
    4805     poly F;
    4806     for(i=1;i<=(size(DD) div 2);i++)
    4807     {
    4808       if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];}
    4809       if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    4810       if(DD[2*i-1]=="movdim"){moverdim=DD[2*i];}
    4811       if(DD[2*i-1]=="version"){version=DD[2*i];}
    4812       if(DD[2*i-1]=="family"){F=DD[2*i];}
    4813     }
    4814     if(F!=0){Fenv=1;}
    4815     list L0;
    4816     if(dimP0(E)==0){L0=2,"Normal";} //  2 es fals pero ha de ser >0 encara que sigui 0
    4817     else
    4818     {
    4819       if(version==0)
    4820       {
    4821           // Computing std(B+E,plex(x,y,x1,..xn)) one can detect if there is a first part
    4822           // independent of parameters giving the variables with dimension 0
    4823        dd=indepparameters(B);
    4824         if (dd==1){d=0; L0=d,string(B);}  // ,determineF(B,F,E)
    4825         else{d=1; L0=2,"Normal";}
    4826       }
    4827       else
    4828       {
    4829         def RH=ringlist(RR);
    4830         def H=RH;
    4831         H[1]=0;
    4832         H[2]=RH[1][2]+RH[2];
    4833         int n=size(H[2]);
    4834         intvec ll;
    4835         for(i=1;i<=n;i++)
    4836         {
    4837           ll[i]=1;
    4838         }
    4839         H[3][1][1]="lp";
    4840         H[3][1][2]=ll;
    4841         def RRH=ring(H);
    4842         setring(RRH);
    4843         ideal BH=imap(RR,B);
    4844         ideal EH=imap(RR,E);
    4845         ideal NH=imap(RR,N);
    4846         if(Fenv==1){poly FH=imap(RR,F);}
    4847         for(i=1;i<=size(EH);i++){BH[size(BH)+1]=EH[i];}
    4848         BH=std(BH);  //  MOLT COSTOS!!!
    4849         ideal G;
    4850         ideal r; poly q;
    4851         for(i=1;i<=size(BH);i++)
    4852         {
    4853           r=factorize(BH[i],1);
    4854           q=1;
    4855           for(j=1;j<=size(r);j++)
    4856           {
    4857             if((pdivi(r[j],NH)[1] != 0) or (equalideals(ideal(NH),ideal(1))))
    4858             {
    4859               q=q*r[j];
    4860             }
    4861           }
    4862           if(q!=1){G[size(G)+1]=q;}
    4863         }
    4864         setring RR;
    4865         def GG=imap(RRH,G);
    4866         ideal GGG;
    4867         if(defined(L0)){kill L0; list L0;}
    4868         for(i=1;i<=size(GG);i++)
    4869         {
    4870           if(indepparameters(GG[i])){GGG[size(GGG)+1]=GG[i];}
    4871         }
    4872         GGG=std(GGG);
    4873         ideal GLM;
    4874         for(i=1;i<=size(GGG);i++)
    4875         {
    4876           GLM[i]=leadmonom(GGG[i]);
    4877         }
    4878         attrib(GLM,"IsSB",1);
    4879         d=dim(std(GLM));
    4880         string antiimi=string(GGG);
    4881         L0=d,antiimi;
    4882         if(d==0)  // d<dim(V(E) \ V(N))?
    4883         {
    4884            " ";string("Anti-image of Special component = ", GGG);
    4885         }
    4886         else
    4887         {
    4888           L0[2]="Normal";
    4889         }
    4890       }
    4891     }
    4892     return(L0);
    4893   }
    48945161
    48955162static proc determineF(ideal A,poly F,ideal E)
     
    49115178  def RRH=ring(H);
    49125179
     5180        //" ";string("Anti-image of Special component = ", GGG);
    49135181
    49145182   setring(RRH);
     
    49195187   ideal M=std(AA+FH);
    49205188   def rh=reduce(EH,M,5);
     5189   //"T_AA="; AA; "T_FH="; FH; "T_EH="; EH; "T_rh="; rh;
    49215190   if(rh==0){env=1;} else{env=0;}
    49225191   setring RR;
    49235192          //L0[3]=env;
     5193    //"T_env="; env;
    49245194    return(env);
    49255195}
    49265196
    49275197
    4928 // locus0(G): Private routine used by locus (the public routine), that
    4929 //               builds the different components.
    4930 // input:     The output G of the grobcov (in generic representation, which is the default option)
    4931 //       Options: The algorithm allows the following options as pair of arguments:
    4932 //               "moverdim", d  : by default moverdim is m=number of variables,
    4933 //                  but it can be set to less values.
    4934 //               "version", v   :  There are two versions of the algorithm. ('version',1) is
    4935 //                  a full algorithm that always distinguishes correctly between 'Normal'
    4936 //                  and 'Special' components, whereas ('version',0) can decalre a component
    4937 //                  as 'Normal' being really 'Special', but is more effective. By default ('version',1)
    4938 //                  is used when the number of variables is less than 4 and 0 if not.
    4939 //                  The user can force to use one or other version, but it is not recommended.
    4940 //               "system", ideal F: if the initial system is passed as an argument.
    4941 //                  This is actually not used.
    4942 //               "comments", c: by default it is 0, but it can be set to 1.
    4943 //                 Usually locus problems have mover coordinates, variables and tracer coordinates.
    4944 //                 The mover coordinates are to be placed as the last variables, and by default,
    4945 //                 its number is dim @P, but as option it can be set to other values.
    4946 // output:
    4947 //         list, the canonical P-representation of the Normal and Non-Normal locus:
    4948 //              The Normal locus has two kind of components: Normal and Special.
    4949 //              The Non-normal locus has two kind of components: Accumulation and Degenerate.
    4950 //              This routine is complemented by locus that calls it in order to eliminate problems
    4951 //              with degenerate points of the mover.
    4952 //         The output components are given as
    4953 //              ((p1,(p11,..p1s_1),type_1,level_1),..,(pk,(pk1,..pks_k),type_k,level_k)
    4954 //         The components are given in canonical P-representation of the subset.
    4955 //              If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level
    4956 //              gives the depth of the component.
    4957 static proc locus0(list GG, list #)
    4958 {
    4959   int te=0;
    4960   int t1=1; int t2=1; int i;
    4961   def RR=basering;
    4962   def Rx=ringlist(RR);
    4963   def P=ring(Rx[1]);
    4964   def RX=Rx;
    4965   RX[1]=0;
    4966   def D=ring(RX);
    4967   def RP=D+P;
    4968   // if(defined(@P)==1){te=1; kill @P; kill @R; kill @RP; }
    4969   // setglobalrings();
    4970   //Options
    4971   list DD=#;
    4972   ideal vmov;
    4973   int dimpar=size(Rx[1][2]);
    4974   int dimvar=size(Rx[2]);
    4975   int nv=nvars(RR);
    4976   int moverdim=nv; //size(ringlist(RR)[1][2]);
    4977    for(i=1;i<=nv;i++)
    4978   {
    4979     //vmov[size(vmov)+1]=var(i+nv-moverdim);
    4980     vmov[size(vmov)+1]=var(i);
    4981   }
    4982   //string("T_moverdim=",moverdim,"; dimvar=",dimvar,"; dimpar=",dimpar,"; vmov=",vmov);
    4983   int version=0;
    4984   //if(nv<4){version=1;}
    4985   int comment=0;
    4986   ideal Fm;
    4987   poly F;
    4988   for(i=1;i<=(size(DD) div 2);i++)
    4989   {
    4990     if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];}
    4991     if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    4992     if(DD[2*i-1]=="version"){version=DD[2*i];}
    4993     if(DD[2*i-1]=="system"){Fm=DD[2*i];}
    4994     if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    4995     if(DD[2*i-1]=="family"){F=DD[2*i];}
    4996   }
    4997   list HHH;
    4998   if (GG[1][1][1]==1 and GG[1][2][1]==1 and GG[1][3][1][1][1]==0 and GG[1][3][1][2][1]==1)
    4999     {return(HHH);}
    5000   list G1; list G2;
    5001   def G=GG;
    5002   list Q1; list Q2;
    5003   int d; int j; int k;
    5004   t1=1;
    5005   for(i=1;i<=size(G);i++)
    5006   {
    5007     attrib(G[i][1],"IsSB",1);
    5008     d=dim(std(G[i][1]));
    5009     if(d==0){G1[size(G1)+1]=G[i];}
    5010     else
    5011     {
    5012       if(d>0){G2[size(G2)+1]=G[i];}
    5013     }
    5014   }
    5015   if(size(G1)==0){t1=0;}  // normal point components
    5016   if(size(G2)==0){t2=0;}  // non-normal point components
    5017   setring(RP);
    5018   if(t1)
    5019   {
    5020     list G1RP=imap(RR,G1);
    5021   }
    5022   else {list G1RP;}
    5023   list P1RP;
    5024   ideal B;
    5025   for(i=1;i<=size(G1RP);i++)
    5026   {
    5027     kill B;
    5028      ideal B;
    5029      for(k=1;k<=size(G1RP[i][3]);k++)
    5030      {
    5031        attrib(G1RP[i][3][k][1],"IsSB",1);
    5032        G1RP[i][3][k][1]=std(G1RP[i][3][k][1]);
    5033        for(j=1;j<=size(G1RP[i][2]);j++)
    5034        {
    5035          B[j]=reduce(G1RP[i][2][j],G1RP[i][3][k][1]);
    5036        }
    5037        B=std(B);
    5038        P1RP[size(P1RP)+1]=list(G1RP[i][3][k][1],G1RP[i][3][k][2],B);
    5039      }
    5040   }  //In P1RP the basis has been reduced wrt the top
    5041    setring(RR);
    5042   ideal h;
    5043   list P1;
    5044   if(t1)
    5045   {
    5046     P1=imap(RP,P1RP);
    5047     for(i=1;i<=size(P1);i++)
    5048     {
    5049       for(j=1;j<=size(P1[i][3]);j++)
    5050       {
    5051         h=factorize(P1[i][3][j],1);
    5052         P1[i][3][j]=h[1];
    5053         for(k=2;k<=size(h);k++)
    5054         {
    5055           P1[i][3][j]=P1[i][3][j]*h[k];
    5056         }
    5057       }
    5058     }
    5059   } //In P1 factors in the basis independent of the parameters have been obtained
    5060   ideal BB; int dd; list NS;
    5061   // "T_AQUI_3";
    5062   for(i=1;i<=size(P1);i++)
    5063   {
    5064      NS=NorSing(P1[i][3],P1[i][1],P1[i][2][1],DD);
    5065      dd=NS[1];
    5066      if(dd==0){P1[i][3]=NS;}  //"Special"; //dd==0
    5067      else{P1[i][3]="Normal";}
    5068   }
    5069   list P2;
    5070   for(i=1;i<=size(G2);i++)
    5071   {
    5072     for(k=1;k<=size(G2[i][3]);k++)
    5073     {
    5074       P2[size(P2)+1]=list(G2[i][3][k][1],G2[i][3][k][2]);
    5075     }
    5076   }
    5077   list l;
    5078   for(i=1;i<=size(P1);i++){Q1[i]=l; Q1[i][1]=P1[i];} P1=Q1;
    5079   for(i=1;i<=size(P2);i++){Q2[i]=l; Q2[i][1]=P2[i];} P2=Q2;
    5080 
    5081   // if(defined(@P)){te=1; kill @P; kill @RP; kill @R;}
    5082   // setglobalrings();
    5083 
    5084   setring(P);
    5085   ideal J;
    5086   if(t1==1)
    5087   {
    5088     def C1=imap(RR,P1);
    5089     def L1=AddLocus(C1);
    5090    }
    5091   else{list C1; list L1; kill P1; list P1;}
    5092   if(t2==1)
    5093   {
    5094     def C2=imap(RR,P2);
    5095     def L2=AddLocus(C2);
    5096   }
    5097   else{list L2; list C2; kill P2; list P2;}
    5098   for(i=1;i<=size(L2);i++)
    5099   {
    5100     J=std(L2[i][2]);
    5101     d=dim(J);
    5102     if(d+1==dimpar)
    5103     {
    5104       L2[i][4]=string("Degenerate",L2[i][4]);
    5105     }
    5106     else{L2[i][4]=string("Accumulation",L2[i][4]);}
    5107   }
    5108 
    5109   list LN;
    5110   if(t1==1)
    5111   {
    5112     for(i=1;i<=size(L1);i++){LN[size(LN)+1]=L1[i];}
    5113   }
    5114   if(t2==1)
    5115   {
    5116     for(i=1;i<=size(L2);i++){LN[size(LN)+1]=L2[i];}
    5117   }
    5118   int tLN=1;
    5119   if(size(LN)==0){tLN=0;}
    5120   setring(RR);
    5121   if(tLN)
    5122   {
    5123    def L=imap(P,LN);
    5124     for(i=1;i<=size(L);i++){if(size(L[i][2])==0){L[i][2]=ideal(1);}}
    5125     //list LL;
    5126     for(i=1;i<=size(L);i++)
    5127     {
    5128       if(typeof(L[i][4])=="list") {L[i][4][1]="Special";}
    5129       l[1]=L[i][2];
    5130       l[2]=L[i][3];
    5131       l[3]=L[i][4];
    5132       l[4]=L[i][5];
    5133       L[i]=l;
    5134     }
    5135   }
    5136   else{list L;}
    5137   // if(te==1){kill @P; kill @RP; kill @R;}
    5138   return(L);
    5139 }
    5140 
    5141 // locus2(G): Private routine used by locus (the public routine), that
    5142 //                builds the different components.
    5143 // input:      The output G of the grobcov (in generic representation, which is the default option)
    5144 //                The ideal F defining the locus problem (G is the grobcov of F and has been
     5198
     5199// locus2(G,F,moverdim,vmov,na):
     5200//                Private routine used by locus (the public routine), that
     5201//                builds the different component, and inputs for locus2
     5202// input:      G= grobcov(S), already computed inside locus
     5203//                F= the ideal defining the locus problem (G is the grobcov of F and has been
     5204//                moverdim=number of mover variables
     5205//                vmov= the ideal of the mover variables
    51455206//                already determined by locus.
    5146 //       Options: The algorithm allows the following options as pair of arguments:
    5147 //                "moverdim", d  : by default moverdim is the number of variables but it can
    5148 //                 be set to minor values.
    5149 //                "version", v   :  There are two versions of the algorithm. ('version',1) is
    5150 //                 a full algorithm that always distinguishes correctly between 'Normal'
    5151 //                 and 'Special' components, whereas ('version',0) can decalre a component
    5152 //                 as 'Normal' being really 'Special', but is more effective. By default ('version',1)
    5153 //                 is used when the number of variables is less than 4 and 0 if not.
    5154 //                 The user can force to use one or other version, but it is not recommended.
    5155 //                 "system", ideal F: if the initial systrem is passed as an argument. This is actually
    5156 //                 not used.
    5157 //                 "comments", c: by default it is 0, but it can be set to 1.
    5158 //                 Usually locus problems have mover coordinates, variables and tracer coordinates.
    5159 //                 The mover coordinates are to be placed as the last variables, and by default,
    5160 //                 its number is equal to the number of parameters of tracer variables, but as option
    5161 //                 it can be set to other values.
     5207//                na= number of parameteres of the locus problem (usually=0);
     5208//                The arguments are determined by locus, and passed to locus2.
    51625209// output:
    51635210//         list, the canonical P-representation of the Normal and Non-Normal locus:
     
    51715218//              If all levels of a class of locus are 1, then the set is locally closed. Otherwise the level
    51725219//              gives the depth of the component.
    5173 // static proc locus2(list G, ideal F, list #)
    5174 static proc locus2(list G, ideal F, list #)
    5175 {
    5176  int st=timer;
    5177  list Snor; list Snonor;
    5178  int d; int i; int j; int mt=0;
    5179  def RR=basering;
    5180  def Rx=ringlist(RR);
    5181  def RP=ring(Rx[1]);
    5182  int tax=1;
    5183  list DD=#;
    5184  list GG=G;
    5185  int n=size(Rx[1][2]);
    5186  int moverdim=n;
    5187  for(i=1;i<=(size(DD) div 2);i++)
    5188   {
    5189      if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i]; mt=1; }
    5190      if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    5191 //    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    5192 //    if(DD[2*i-1]=="version"){version=DD[2*i];}
    5193 //    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    5194 //    if(DD[2*i-1]=="grobcov"){GG=DD[2*i];}
    5195 //    if(DD[2*i-1]=="anti-image"){tax=DD[2*i];}
    5196   }
    5197   if(mt==0){DD[size(DD)+1]="moverdim"; DD[size(DD)+1]=moverdim;}
    5198  for(i=1;i<=size(GG);i++)
    5199  {
    5200    attrib(GG[i][1],"IsSB",1);
    5201    GG[i][1]=std(GG[i][1]);
    5202    d=dim(GG[i][1]);
    5203    if(d==0)
     5220static proc locus2(list G, ideal F, int moverdim, ideal vmov, int na)
     5221{
     5222   int st=timer;
     5223   list Snor; list Snonor;
     5224   int d; int i; int j; //int mt=0;
     5225   def RR=basering;
     5226   def Rx=ringlist(RR);
     5227   def RP=ring(Rx[1]);
     5228   def LP=ringlist(RP);
     5229   int nax=size(LP[2]);
     5230   int nx=nax-na;
     5231   int nv=moverdim;
     5232   int tax=1;
     5233   list GG=G;
     5234   int n=size(Rx[1][2]);
     5235   for(i=1;i<=size(GG);i++)
    52045236   {
    5205      for(j=1;j<=size(GG[i][3]);j++)
    5206      {
    5207        Snor[size(Snor)+1]=GG[i][3][j];
    5208      }
    5209    }
    5210    else
    5211    {
    5212      if(d>0)
     5237     attrib(GG[i][1],"IsSB",1);
     5238     GG[i][1]=std(GG[i][1]);
     5239     d=dim(GG[i][1]);
     5240     if(d==0)
    52135241     {
    52145242       for(j=1;j<=size(GG[i][3]);j++)
    52155243       {
    5216          Snonor[size(Snonor)+1]=GG[i][3][j];
    5217       }
     5244         Snor[size(Snor)+1]=GG[i][3][j];
     5245       }
     5246     }
     5247     else
     5248     {
     5249       if(d>0)
     5250       {
     5251         for(j=1;j<=size(GG[i][3]);j++)
     5252         {
     5253           Snonor[size(Snonor)+1]=GG[i][3][j];
     5254        }
     5255       }
    52185256     }
    52195257   }
    5220  }
    5221  int tnor=size(Snor); int tnonor=size(Snonor);
    5222  setring RP;
    5223  list SnorP;
    5224  list SnonorP;
    5225  if(tnor)
    5226  {
    5227    SnorP=imap(RR,Snor);
    5228    st=timer;
    5229    SnorP=LCUnionN(SnorP);
    5230  }
    5231  if(tnonor)
    5232  {
    5233    SnonorP=imap(RR,Snonor);
    5234    SnonorP=LCUnionN(SnonorP);
    5235  }
    5236  setring RR;
    5237  ideal C;  list N;  list BAC; list AI;
    5238  list NSC; list DAC;
    5239  list L;
    5240  ideal B;
    5241  int k;
    5242  int j0; int k0; int te;
    5243  poly kkk=1;
    5244  ideal AI0;
    5245  int dimP;
    5246 
    5247  if(tnor)
    5248  {
    5249    Snor=imap(RP,SnorP);
    5250    for(i=1;i<=size(Snor);i++)
     5258   //"T_Snor="; Snor;
     5259   //"T_Snonor="; Snonor;
     5260   int tnor=size(Snor); int tnonor=size(Snonor);
     5261   setring RP;
     5262   list SnorP;
     5263   list SnonorP;
     5264   if(tnor)
    52515265   {
    5252      C=Snor[i][1];
    5253      N=Snor[i][2];
    5254      dimP=DimPar(C);
    5255      AI=NS(F,G,C,N,DD);
    5256      AI0=AI[2];
    5257      if(size(AI)==3){AI[2]="normal"; AI[3]=ideal(0);}
    5258      else
     5266     SnorP=imap(RR,Snor);
     5267     st=timer;
     5268     SnorP=LCUnionN(SnorP);
     5269   }
     5270   if(tnonor)
     5271   {
     5272     SnonorP=imap(RR,Snonor);
     5273     SnonorP=LCUnionN(SnonorP);
     5274   }
     5275   //"T_SnorP after LCUnion="; SnorP;
     5276   // "T_SnonorP  after LCUnion="; SnonorP;
     5277   setring RR;
     5278   ideal C;  list N;  list BAC; list AI;
     5279   list NSC; list DAC;
     5280   list L;
     5281   ideal B;
     5282   int k;
     5283   int j0; int k0; int te;
     5284   poly kkk=1;
     5285   ideal AI0;
     5286   int dimP;
     5287
     5288   if(tnor)
     5289   {
     5290     Snor=imap(RP,SnorP);
     5291     for(i=1;i<=size(Snor);i++)
    52595292     {
    5260        if(AI[1]>=dimP)
    5261        {
    5262          AI[2]="Normal";
    5263          if(tax>0){AI[3]=AI0;}
    5264        }
    5265        else
    5266        {
    5267            if(AI[1]<n-1){AI[2]="Special"; AI[3]=AI0;}
    5268        }
     5293       C=Snor[i][1];
     5294       N=Snor[i][2];
     5295       dimP=DimPar(C,nax,nx);
     5296        //"T_G="; G;
     5297       AI=NS(F,G,C,N,moverdim,na,vmov,dimP);
     5298       Snor[i][size(Snor[i])+1]=AI;
    52695299     }
    5270      Snor[i][size(Snor[i])+1]=AI;
     5300     for(i=1;i<=size(Snor);i++)
     5301     {
     5302       L[size(L)+1]=Snor[i];
     5303     }
     5304    }
     5305    ideal AINN;
     5306   if(tnonor)
     5307   {
     5308     Snonor=imap(RP,SnonorP);
     5309     //"T_Snonor="; Snonor;
     5310     //"T_G="; G;
     5311     for(i=1;i<=size(Snonor);i++)
     5312     {
     5313       DAC=DimComp(Snonor[i][1],nax,nx);
     5314       Snonor[i][size(Snonor[i])+1]=DAC;
     5315     }
     5316     for(i=1;i<=size(Snonor);i++)
     5317     {
     5318       L[size(L)+1]=Snonor[i];
     5319     }
    52715320   }
    5272    for(i=1;i<=size(Snor);i++)
    5273    {
    5274      L[size(L)+1]=Snor[i];
    5275    }
    5276   }
    5277   ideal AINN;
    5278  if(tnonor)
    5279  {
    5280    Snonor=imap(RP,SnonorP);
    5281    for(i=1;i<=size(Snonor);i++)
    5282    {
    5283      DAC=dimComp(Snonor[i][1]);
    5284      Snonor[i][size(Snonor[i])+1]=DAC;
    5285 //      if(tax==1)
    5286 //      {
    5287 //        AINN=AISnonor(Snonor[i][1],GG,DD);
    5288 //        Snonor[i][3][3]=AINN;
    5289 //       // AQUI ES ON AFEGIM LA ANTIMATGE DE LES NONOR
    5290 //      }
    5291    }
    5292    for(i=1;i<=size(Snonor);i++)
    5293    {
    5294      L[size(L)+1]=Snonor[i];
    5295    }
    5296  }
    5297  return(L);
    5298 }
    5299 
    5300 // Auxiliary algorithm of locus2.
     5321  return(L);
     5322}
     5323
     5324// Auxilliary algorithm of locus2.
    53015325//           The algorithm searches the basis corresponding to C, in the grobcov.
    53025326//           It reduces the basis modulo the component.
     
    53045328//           For each hole of the component
    53055329//              it searches the segment where the hole is included
    5306 //              and select form its basis the polynomials
     5330//              and selects the polynomials from its basis
    53075331//              only dependent on the variables.
    53085332//              These polynomials are non-null in an open set of
     
    53205344//           taxonomy is not determined, and (n,'normal', 0) is returned as third
    53215345//           element of the component. (n is the dimension of the space).
    5322 static proc NS(ideal F,list G, ideal C, list N, list #)
     5346static proc NS(ideal F,list G, ideal C, list N, int nv, int na,ideal vmov,int dimC)
    53235347{
    53245348  // Initializing and defining rings
    5325   int i; int j; int k; int te; int j0;int k0; int m;
    5326   list DD=#;
    5327   def RR=basering;
     5349   int i; int j; int k; int te; int j0;int k0; int m;
     5350   def RR=basering;
     5351   def Lax_uv=ringlist(RR);
     5352   Lax_uv[3][1][1]="lp";
     5353   int nax=size(Lax_uv[1][2]);
     5354   int nuv=size(Lax_uv[2]);
     5355   int nx=nax-na;
     5356   int nu=nuv-nv;
     5357   //"Lax_uv="; Lax_uv;
     5358    def Lax=Lax_uv[1];
     5359    def Qax=ring(Lax);                                  // ring Q[a,x]
     5360   //"T_Lax="; Lax;
     5361   def La_x=LaxToLa_x(Lax,nx);                   // ring Q[a][x]
     5362   //"T_La_x="; La_x;
     5363   def  La_v=Lax_uvToLa_v(Lax_uv,na,nv);    // ring Q[a][v]
     5364   //"T_La_v="; La_v;
     5365   def  Lxuav=Lax_uvToLxuav(Lax_uv,na,nv);
     5366  //"T_Lxuav="; Lxuav;
     5367
     5368
     5369  // old rings
    53285370  def Rx=ringlist(RR);
    53295371  def Lx=Rx;
    5330   def P=ring(Rx[1]);               // ring Q[bx]
     5372  def P=ring(Rx[1]);                 // ring Q[a,x]]
    53315373  Lx[1]=0;
    5332   def D=ring(Lx);                   // ring Q[bu]
    5333   def PR0=P+D;
     5374  def D=ring(Lx);                     // ring Q[u,v]
     5375  def PR0=P+D;                       // ring Q[a,x,u,v]
    53345376  def PRx=ringlist(PR0);
    53355377  PRx[3][2][1]="lp";
    5336   def PR=ring(PRx);               // ring Q[bx,bu]  in lex order
    5337   int n=size(Rx[1][2]);           // number of parameters
    5338   def nv=nvars(RR);               // number of variables
    5339   int moverdim=n;                // number of mover variables
    5340   ideal vmov;                         // mover variables   (bw)
    5341   for(i=1;i<=(size(DD) div 2);i++)
    5342   {
    5343     if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    5344     if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];}
    5345 //    if(DD[2*i-1]=="version"){version=DD[2*i];}
    5346 //    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    5347 //    if(DD[2*i-1]=="grobcov"){GG=DD[2*i];}
    5348 //    if(DD[2*i-1]=="anti-image"){tax=DD[2*i];}
    5349   }
    5350 
    5351   //string("T_moverdim= ", moverdim);
    5352   // vmov = the last vmov variables = mover variables
    5353   for(i=1;i<=moverdim;i++)
    5354   {
    5355     vmov[size(vmov)+1]=var(i+nv-moverdim);
    5356   }
    5357   //string("T_nv=",nv,"  moverdim=",moverdim);
    5358 //  for(i=1;i<=nv;i++)
    5359 //  {
    5360 //    vmov[size(vmov)+1]=var(i);
    5361 //  }
    5362 
     5378  // "T_PRx="; PRx;
     5379  def PR=ring(PRx);                 // ring Q[a,x,u,v]  in lex order
     5380  // end of old rings
     5381
     5382  for(i=1;i<=nv;i++)
     5383  {
     5384    vmov[size(vmov)+1]=var(i+nuv-nv);
     5385  }
     5386  //string("T_nv=",nv,"  moverdim=",nv);
    53635387  int ddeg; int dp;
    53645388  list LK;
     
    53665390  for(i=1;i<=nv;i++){bu[i]=var(i);}
    53675391  ideal mv;
    5368 // ULL
    5369 // ideal of mover varibles
    5370 //  for(i=1;i<=nv;i++){mv[size(mv)+1]=var(i);}
    53715392   for(i=1;i<=nv;i++){mv[size(mv)+1]=var(i);}
    53725393
     
    53855406   if(te==1){"ERROR";}
    53865407   def B=G[j0][2];  // Aixo aniria be per les nonor
    5387 
    5388    // Searching the elements in Q[v_m]  on basis different from B that are nul there
     5408   //"T_B=G[j0][2]="; B;
     5409   //string("T_k0=",k0," G[",j0,"]="); G[j0];
     5410
     5411   // Searching the elements in Q[v_m]  on basis differents from B that are nul there
    53895412   // and cannot become 0 on the antiimage of B. They are placed on NoNul
    53905413   list NoNul;
     
    54405463  }
    54415464
    5442    // Adding F to B
    5443    for(i=1;i<=size(F);i++)
    5444    {
    5445      B[size(B)+1]=F[i];
    5446    }
     5465  // Adding F to B
     5466  for(i=1;i<=size(F);i++)
     5467  {
     5468    B[size(B)+1]=F[i];
     5469  }
    54475470
    54485471  def E=NoNul;
     
    54575480  setring(RR);
    54585481  // if(n+nv>10 or ddeg>=10){LK=n,ideal(0),"normal"; return(LK);}
    5459   if(ddeg>=10){LK=n,ideal(0),"normal"; return(LK);}  // 8 instead of 10 ?
     5482  if(ddeg>=10){LK=nv,ideal(0),"normal"; return(LK);}  // 8 instead of 10 ?
    54605483
    54615484  // Reducing  basis B modulo C  in the ring PR      lex(x,u)
    5462   setring(PR);
    5463   def buPR=imap(RR,bu);
    5464   def EE=imap(RR,E);
    5465   ideal BR;
    5466   def BB=imap(RR,B);
     5485  // setring(PR);
     5486
     5487  def Qxuav=ring(Lxuav);
     5488  setring(Qxuav);
     5489  def BR=imap(RR,B);
     5490  ideal vamov;
     5491  for(i=nx+nu+1;i<=nax+nuv;i++){vamov[size(vamov)+1]=var(i);}
     5492
     5493  //BR=std(BR);
    54675494  def CC=imap(RR,C);
     5495  for(i=1;i<=size(CC);i++){BR[size(BR)+1]=CC[i];}
     5496  BR=std(BR);
     5497 // for(i=1;i<=size(CC);i++){BR[size(BR)+1]=CC[i];}
    54685498  attrib(CC,"IsSB",1);
    5469   CC=std(CC);
    5470   for(i=1;i<=size(BB);i++)
    5471   {
    5472     BR[i]=reduce(BB[i],CC);
    5473   }
    5474 
    5475   // Computing the GB(B) in the ring PR     lex(x,u)
    5476   BR=std(BR);
    5477   ideal BRP;
    5478 
    5479  // Eliminating the polynomials in B containing variables v_t
    5480   for(i=1;i<=size(BR);i++)
    5481   {
    5482     if(subset(variables(BR[i]),buPR )){BRP[size(BRP)+1]=BR[i];}
    5483   }
    5484   BR=BRP;
    5485 
    5486   // Factoring and eliminating the factors in N
    5487   list H;
    5488   ideal BP;
    5489   poly f;
    5490   int tj;
    5491   //ideal vv=imap(RR,bu);
    5492   for(i=1;i<=size(BR);i++)
    5493   {
    5494     f=1;
    5495     H=factorize(BR[i]);
    5496     if((subset(variables(H[1][2]),buPR)) and (size(H[1]))==2){BP[size(BP)+1]=H[1][2];}
    5497     else
    5498     {
    5499       for(j=1;j<=size(H[1]);j++)
    5500       {
    5501         tj=subset(H[1][j],EE);
    5502         if((subset(variables(H[1][j]),buPR)) and (tj==0)){f=f*H[1][j];}
    5503       }
    5504       if(leadcoef(f)!=f){BP[size(BP)+1]=f;}
    5505    }
    5506   }
    5507 
    5508   // Determining the GB basis of B in Q[u]
    5509   BP=std(BP);
    5510   // Determining the dimension of the anti-image
    5511   //   in the ring D
    5512   setring D;
    5513   def KK=imap(PR,BP);
    5514   KK=std(KK);
    5515   ideal KKM;
    5516   ideal vmovM=imap(RR,vmov);
    5517 
    5518   // selecting the element in Q[w]
    5519   for(i=1;i<=size(KK);i++)
    5520   {
    5521     if(subset(variables(KK[i]),vmovM))
    5522     {
    5523       KKM[size(KKM)+1]=KK[i];
    5524     }
    5525   }
    5526 
    5527   // Determining the dimensions
    5528   list AIM=dimM(KKM,nv,moverdim);
    5529   //int d2=AIM[1];
    5530   def AAIM=AIM[2];
     5499  ideal AIM;
     5500 // "T_BR="; BR;
     5501  for(i=1;i<=size(BR);i++){if(subset(variables(BR[i]),vamov)){AIM[size(AIM)+1]=BR[i];}}
     5502  //"T_AIM="; AIM;
     5503
     5504  list La_v0=imap(RR,La_v);
     5505  def Qa_v=ring(La_v0);
     5506  setring Qa_v;
     5507  def AIMa_v=imap(Qxuav,AIM);
     5508  AIMa_v=std(AIMa_v);
     5509  int dimAIM=dim(AIMa_v);
     5510  //"T_AIMa_v="; AIMa_v;
     5511  //string("T_dimAIM=",dimAIM);
    55315512  setring(RR);
    5532   //def BBR=imap(D,KK);
    5533   def LKK=imap(D,AIM);
    5534   //LKK=d2,string(BBR);
    5535   return(LKK);
    5536 }
    5537 
    5538 // Auxiliary algorithm of locus2.
    5539 //           The algorithm searches the basis corresponding to C, in the grobcov.
    5540 //           It reduces the basis modulo the component.
    5541 //           The result is the reduced basis BR.
    5542 //           for each hole of the component
    5543 //              it searches the segment where the hole is included
    5544 //              and select form its basis the polynomials
    5545 //              only dependent on the variables.
    5546 //              These polynomials are non-null in an open set of
    5547 //              the component, and are included in the list NoNul of non-null factors
    5548 // input: C: The top of the the segment of a non-normal locus
    5549 //           G the grobcov of F
    5550 // output: B: the anti-image of the input segment
    5551 static proc AISnonor(ideal C,list G, list #)
    5552 {
    5553   // Initializing and defining rings
    5554   int i; int j; int k; int te; int j0;int k0; int m;
    5555   list DD=#;
    5556   def RR=basering;
    5557   def Rx=ringlist(RR);
    5558   def Lx=Rx;
    5559   def P=ring(Rx[1]);               // ring Q[bx]
    5560   Lx[1]=0;
    5561   def D=ring(Lx);                   // ring Q[bu]
    5562   def PR0=P+D;
    5563   def PRx=ringlist(PR0);
    5564   PRx[3][2][1]="lp";
    5565   def PR=ring(PRx);               // ring Q[bx,bu]  in lex order
    5566   int n=size(Rx[1][2]);           // number of parameters
    5567   def nv=nvars(RR);               // number of variables
    5568   int moverdim=n;                // number of parameters
    5569   ideal vmov;                         // mover variables   (bw)
    5570 
    5571   for(i=1;i<=(size(DD) div 2);i++)
    5572   {
    5573     if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    5574     if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];}
    5575 //    if(DD[2*i-1]=="version"){version=DD[2*i];}
    5576 //    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    5577 //    if(DD[2*i-1]=="grobcov"){GG=DD[2*i];}
    5578 //    if(DD[2*i-1]=="anti-image"){tax=DD[2*i];}
    5579   }
    5580   for(i=1;i<=moverdim;i++)
    5581   {
    5582     vmov[size(vmov)+1]=var(i+nv-moverdim);
    5583   }
    5584 
    5585   int ddeg; int dp;
    5586   list LK;
    5587   ideal bu;                               // ideal of all variables   (u)
    5588   for(i=1;i<=nv;i++){bu[i]=var(i);}
    5589   ideal mv;                            // ideal of mover varibles
    5590   for(i=1;i<=nv;i++){mv[size(mv)+1]=var(i);}
    5591 
    5592    // Searching the basis associated to Nonor
    5593    j=2; te=1;
    5594    while((te) and (j<=size(G)))
    5595    {
    5596       k=1;
    5597       while((te) and (k<=size(G[j][3])))
    5598       {
    5599         if (equalideals(C,G[j][3][k][1])){j0=j; k0=k; te=0;}
    5600         k++;
    5601       }
    5602       j++;
    5603    }
    5604    if(te==1){"ERROR";}
    5605    def B=G[j0][2];
    5606    return(B);
    5607 }
    5608 
    5609 // NOT USED
    5610 // Auxiliary algorithm of locus2.
    5611 //           The algorithm searches the basis corresponding to C, in the grobcov.
    5612 //           It reduces the basis modulo the component.
    5613 //           The result is the reduced basis BR.
    5614 //           for each hole of the component
    5615 //              it searches the segment where the hole is included
    5616 //              and select form its basis the polynomials
    5617 //              only dependent on the variables.
    5618 //              These polynomials are non-null in an open set of
    5619 //              the component, and are included in the list NoNul of non-null factors
    5620 // input: F: the ideal of the locus problem
    5621 //           G the grobcov of F
    5622 //           C the top of a component of normal points
    5623 //           N the holes of the component
    5624 // output: (d,tax,a)
    5625 //           where d is the dimension of the anti-image
    5626 //           a is the anti-image of the component and
    5627 //           tax is the taxonomy \"Normal\" if d is equal to the dimension of C
    5628 //           and \"Special\" if it is smaller.
    5629 static proc NS5(ideal F,list G, ideal C, list N)
    5630 {
    5631   // Initializing and defining rings
    5632   int i; int j; int k; int te; int j0;int k0; int m;
    5633   def RR=basering;
    5634   def Rx=ringlist(RR);
    5635   def Lx=Rx;
    5636   def P=ring(Rx[1]);
    5637   Lx[1]=0;
    5638   def D=ring(Lx);
    5639   def RP=D+P;
    5640   def PR0=P+D;
    5641   def PRx=ringlist(PR0);
    5642   PRx[3][2][1]="lp";
    5643   def PR=ring(PRx);
    5644   int n=size(Rx[1][2]);                                     // dimension of the tracer space
    5645   def nv=nvars(RR);                                        // number of variables
    5646   ideal v; for(i=1;i<=nv;i++){v[i]=var(i);}        // variables
    5647   int nm;                                                        // number of mover variables
    5648   ideal vmov;                                                  // mover variables
    5649   if(nv>=n){nm=n;}
    5650   else{nm=nv;}
    5651   for(i=1;i<=nm;i++)
    5652   {
    5653     vmov[size(vmov)+1]=var(i+nv-nm);
    5654   }
    5655 
    5656   // Searching the basis associated to C
    5657    j=2; te=1;
    5658    while((te) and (j<=size(G)))
    5659    {
    5660       k=1;
    5661       while((te) and (k<=size(G[j][3])))
    5662       {
    5663         if (equalideals(C,G[j][3][k][1])){j0=j; k0=k; te=0;}
    5664         k++;
    5665       }
    5666       j++;
    5667    }
    5668    if(te==1){"ERROR";}
    5669    def B=G[j0][2];     // basis associated to C
    5670 
    5671   // Searching the elements in Q[vmov]  on bases different from B that are nul there
    5672   // and cannot become 0 on the antiimage of C. They are placed on NoNul
    5673   list NoNul;          // elements in Q[vmov]  on bases different from B
    5674   ideal BNoNul;
    5675   ideal covertop;   // basis of the segment whose top is a hole of our segment
    5676   int te1;
    5677   for(i=1;i<=size(N);i++)
    5678   {
    5679     j=2; te=1;
    5680     while(te and j<=size(G))
    5681     {
    5682       if(j!=j0)
    5683       {
    5684         k=1;
    5685         while(te and k<=size(G[j][3]))
    5686         {
    5687           covertop=G[j][3][k][1];
    5688           if(equalideals(covertop,N[i]))
    5689           {
    5690             te=0; te1=1; BNoNul=G[j][2];
    5691           }
    5692           else
    5693           {
    5694             if(redPbasis(covertop,N[i]))
    5695             {
    5696               te=0; te1=1; m=1;
    5697               while( te1 and m<=size(G[j][3][k][2]) )
    5698               {
    5699                 if(equalideals(G[j][3][k][2][m] ,N[i] )==1){te1=0;}
    5700                 m++;
    5701               }
    5702             }
    5703             if(te1==1){ BNoNul=G[j][2];}
    5704           }
    5705           k++;
    5706         }
    5707 
    5708         if((te==0) and (te1==1))
    5709         {
    5710           // Selecting the elements independent of the parameters,
    5711           // They will be non null on the segment
    5712           for(m=1;m<=size(BNoNul);m++)
    5713           {
    5714             if(indepparameterspoly(BNoNul[m]))
    5715             {
    5716                NoNul[size(NoNul)+1]=BNoNul[m];
    5717             }
    5718          }
    5719        }
    5720      }
    5721      j++;
    5722    }
    5723  }
    5724   def NN=NoNul;
    5725   poly kkk=1;
    5726   if(size(NN)==0){NN[1]=kkk;}
    5727 
    5728  // Union of F and B
    5729  for(i=1;i<=size(F);i++)
    5730  {
    5731    B[size(B)+1]=F[i];
    5732  }
    5733 
    5734   // Reducing  basis B modulo C  in the ring PR
    5735   setring(PR);
    5736   def vv=imap(RR,v);
    5737   ideal BBR;
    5738   def BB=imap(RR,B);
    5739   def CC=imap(RR,C);
    5740   attrib(CC,"IsSB",1);
    5741   CC=std(CC);
    5742   for(i=1;i<=size(BB);i++)
    5743   {
    5744     BBR[i]=reduce(BB[i],CC);
    5745   }
    5746 
    5747   // Selecting the elements of the ideal in Q[v]
    5748   // setring RP;
    5749   ideal Bv;
    5750   // def Bv=imap(PR,BBR);
    5751   //ideal vvv=imap(RR,v);
    5752   // Bv=std(Bv);
    5753   ideal BR;
    5754   for(i=1;i<=size(BBR);i++)
    5755   {
    5756     if(subset(variables(BBR[i]),vv)){BR[size(BR)+1]=BBR[i];}
    5757   }
    5758   // setring PR;
    5759   // def BR=imap(RP,BRv);
    5760   def NNN=imap(RR,NN);
    5761 
    5762   // Factoring and eliminating the factors in N
    5763   list H;
    5764   ideal BP;
    5765   poly f;
    5766   int tj;
    5767   // ideal vv=imap(RR,v);
    5768   for(i=1;i<=size(BR);i++)
    5769   {
    5770     f=1;
    5771     H=factorize(BR[i]);
    5772     if((subset(variables(H[1][2]),vv)) and (size(H[1]))==2){BP[size(BP)+1]=H[1][2];}
    5773     else
    5774     {
    5775       for(j=1;j<=size(H[1]);j++)
    5776       {
    5777         tj=subset(H[1][j],NNN);
    5778         if((subset(variables(H[1][j]),vv)) and (tj==0)){f=f*H[1][j];}
    5779       }
    5780       if(leadcoef(f)!=f){BP[size(BP)+1]=f;}
    5781    }
    5782   }
    5783 
    5784   //  Obtaining Am, the anti-image of C on vmov
    5785   setring D;
    5786   ideal A=imap(PR,BP);
    5787   ideal vm=imap(RR,vmov);
    5788   ideal Am;
    5789   for(i=1;i<=size(A);i++)
    5790   {
    5791     if (subset(variables(A[i]),vm)){Am[size(Am)+1]=A[i];}
    5792   }
    5793 
    5794   list AIM=dimM(Am,nv,nm);
    5795   setring(RR);
    5796   def LAM=imap(D,AIM);
    5797   return(LAM);
    5798 }
    5799 
    5800 static proc dimM(ideal KKM, int nv, int nm)
     5513  def AIMRR=imap(Qa_v,AIMa_v);
     5514  string TaxComp;
     5515  if(dimAIM==dimC){TaxComp="Normal";}
     5516  else{TaxComp="Special";}
     5517  list NSA=dimAIM,TaxComp,AIMRR;
     5518  return(NSA);
     5519}
     5520
     5521static proc DimM(ideal KKM, int na, int nv)
    58015522{
    58025523  def RR=basering;
     
    58045525  int i;
    58055526  def Rx=ringlist(RR);
    5806   for(i=1;i<=nm;i++)
     5527
     5528  for(i=1;i<=nv;i++)
    58075529  {
    58085530    L[i]=Rx[2][nv-nm+i];
     
    58135535  Rx[3][1][2]=iv;
    58145536   def DM=ring(Rx);
     5537  //"Rx="; Rx;
    58155538  setring(DM);
    58165539  ideal KKMD=imap(RR,KKM);
     
    58215544  def KAIM=imap(DM,KKMD);
    58225545  list LAIM=d,KAIM;
     5546 // "T_LAIM="; LAIM;
    58235547  return(LAIM);
    58245548}
     
    58585582  ideal vpar;
    58595583  ideal vvar;
     5584  //"T_n="; n;
     5585   //"T_nv="; nv;
    58605586  for(i=1;i<=n;i++){vpar[size(vpar)+1]=par(i);}
    58615587  for(i=1;i<=nv;i++){vvar[size(vvar)+1]=var(i);}
     5588  //string("T_vpar = ", vpar," vvar = ",vvar);
    58625589  def P=ring(Rx[1]);
    58635590  Rx[1]=0;
     
    58705597  //string("T_vvpar = ",vvpar);
    58715598  ideal B=std(FF);
     5599  //"T_B="; B;
    58725600  ideal Bel;
     5601  //"T_vvpar="; vvpar;
    58735602  for(i=1;i<=size(B);i++)
    58745603  {
    58755604    if(subset(variables(B[i]),vvpar)) {Bel[size(Bel)+1]=B[i];}
    58765605  }
     5606  //"T_Bel="; Bel;
    58775607  list H;
    58785608  list FH;
     
    59045634//  locus(F):  Special routine for determining the locus of points
    59055635//                 of  geometrical constructions.
    5906 //  input:      The ideal of the locus equations
     5636//  input:      The ideal of the locus equations defined in the
     5637//                 ring Q[a1,..,ap,x1,..xn][u1,..um,v1,..vn]
    59075638//  output:
    59085639//          The output components are given as
     
    59115642//               The third element tax is:
    59125643//                 for normal point components, tax=(d,taxonomy,anti-image)
    5913 //                    being d=dimension of the anti-image on the mover varaibles,
     5644//                    being d=dimension of the anti-image on the mover variables,
    59145645//                           taxonomy='Normal'  or  'Special', and
    59155646//                           anti-image=ideal of the anti-image over the mover variables
    5916 //                                which are taken to be the last n variables.
     5647//                                which by default are taken to be the last n variables.
    59175648//                 for non-normal point components, tax =(d,taxonomy)
    59185649//                    being d=dimension of the component  and
     
    59225653//            Normal component:
    59235654//            - each point in the component has 0-dimensional anti-image.
    5924 //            - the anti-image in the mover coordinates is equal to the dimension of the component
     5655//            - the anti-image in the mover coordinates is equal to the dimension of the component.
    59255656//            Special component:
    59265657//            - each point in the component has 0-dimensional anti-image.
    5927 //            - the anti-image in the mover coordinates is smaller than the dimension of the component
     5658//            - the anti-image on the mover variables is smaller than the dimension of the component.
    59285659//          The non-normal locus has two kind of components: Accumulation and Degenerate.
    59295660//           Accumulation points:
     
    59365667//           taxonomy is not determined, and (n,'normal', 0) is returned as third
    59375668//           element of the component. (n is the dimension of the space).
    5938 
    59395669proc locus(ideal F, list #)
    5940 "USAGE:  locus(ideal F[,options])
     5670"USAGE:  locus(ideal F [,options])
    59415671        Special routine for determining the locus of points of
    59425672        a geometrical construction.
    5943 INPUT:  The input ideal must be the set equations
    5944         defining the locus. Calling sequence:
    5945         locus(F[,options]);
    5946         The input ring must be a parametrical ideal
    5947         in Q[x1,..,xn][u1,..,um],
    5948         (x=tracer variables, u=mover variables).
    5949         The tracer variables are x1,..xn = dim(x-space).
    5950         By default, the mover variables are the last n u-variables.
     5673INPUT:  The input ideal must be the ideal of the set equations
     5674        defining the locus, defined in the ring
     5675        ring Q(0,a1,..,ap,x1,..xn)(u1,..um,v1,..vn),lp;
     5676        Calling sequence:
     5677        locus(F [,options]);
     5678        a=fixed parameters,x=tracer variables, u=auxiliary variables, v=mover variables.
     5679        The parameters a are optative. If they are used, then the option \"numpar\=,np
     5680        must be declared, being np the number of fixed parameters.
     5681        The tracer variables are x1,..xn, where n is the dimension of the space.
     5682        By default, the mover variables are the last n variables.
    59515683        Its number can be forced by the user to the last
    59525684        k variables by adding the option \"moverdim\",k.
    59535685        Nevertheless, this option is recommended only
    59545686        to experiment, and can provide incorrect taxonomies.
     5687        The remaining variables are auxiliary.
    59555688OPTIONS: An option is a pair of arguments: string, integer.
    59565689        To modify the default options, pairs of arguments
     
    59595692        pair of arguments:
    59605693
    5961         \"moverdim\", k  to restrict the mover-variables of the
    5962         anti-image to the last k variables.
    5963         By defaulat k is equal to the last n u-variables,
     5694        \"numpar\", np  in order to consider the first np parameters of the ring
     5695        to be fixed parameters of the locus, being the tracer variables
     5696        the remaining parameters.
     5697        To be used for a paramteric locus. (New in release N12).
     5698
     5699        \"moverdim\", k  to force the mover-variables to be the last
     5700         k variables. This determines the antiimage and its dimension.
     5701        By defaulat k is equal to the last n variables,
    59645702        We can experiment with a different value,
    59655703        but this can produce an error in the character
    59665704         \"Normal\" or \"Special\" of a locus component.
    59675705
    5968         \"grobcov\", L, where L is the list of a previous computed grobcov(F).
     5706        \"grobcov\", G, where G is the list of a previous computed grobcov(F).
    59695707        It is to be used when we modify externally the grobcov,
    59705708        for example to obtain the real grobcov.
    59715709
    5972         \"comments\", c: by default it is 0, but it can be set
    5973         to 1.
     5710        \"comments\", c: by default it is 0, but it can be set to 1.
    59745711RETURN: The output is a list of the components:
    59755712        ((p1,(p11,..p1s_1),tax_1), .., (pk,(pk1,..pks_k),tax_k)
     
    59985735           - each point in the component has 0-dimensional
    59995736              anti-image.
    6000            - the anti-image in the mover coordinates is smaller
    6001               than the dimension of the component
     5737           - the anti-image in the mover coordinates has dimension
     5738              smaller than the dimension of the component
    60025739        The non-normal locus has two kind of components:
    60035740          Accumulation and Degenerate.
     
    60135750         then the taxonomy is not determined, and (n,'normal', 0)
    60145751         is returned as third element of the component. (n is the
    6015          dimension of the space).
     5752         dimension of the tracer space).
    60165753
    60175754       Given a parametric ideal F representing the system F
     
    60315768EXAMPLE: locus; shows an example"
    60325769{
    6033   int tes=0; int i;  int m; int mm; int n;
     5770  int tes=0; int i;  int m; int mm; // int n;
    60345771  def RR=basering;
    60355772  list GG;
    60365773  //Options
    60375774  list DD=#;
    6038   ideal vmov;
    6039   int nv=nvars(RR);                                  // number of variables
    6040  // int moverdim=nv; //size(ringlist(RR)[1][2]);   // number of parameters
    6041   int moverdim=size(ringlist(RR)[1][2]);   // number of parameters
    6042     if(moverdim>nv){moverdim=nv;}
    6043   // for(i=1;i<=nv;i++){vmov[size(vmov)+1]=var(i);}
    6044   int version=2;
    6045   // if(nv<4){version=2;}
     5775  int nax=size(ringlist(RR)[1][2]);              // number of parameters + tracer variables
     5776  int nuv=nvars(RR);                                 // number of variables
     5777  int na=0; int nx=nax;
     5778  int moverdim=nx;                                  // number of tracer variables
     5779  if(moverdim>nuv){moverdim=nuv;}
     5780//  int version=2;
    60465781  int comment=0;
    60475782  int tax=1;
     
    60495784  for(i=1;i<=(size(DD) div 2);i++)
    60505785  {
    6051     if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];}
    6052     if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    6053     // if(DD[2*i-1]=="version"){version=DD[2*i];}
     5786    if(DD[2*i-1]=="numpar"){na=DD[2*i]; nx=nax-na; moverdim=nx;}
    60545787    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    60555788    if(DD[2*i-1]=="grobcov"){GG=DD[2*i];}
    60565789  }
    6057   // if(version != 2){string("versions 0 and 1 are not available"); version=2;}
    6058   for(i=1;i<=moverdim;i++){vmov[size(vmov)+1]=var(i+nv-moverdim);}
     5790  for(i=1;i<=(size(DD) div 2);i++)
     5791  {
     5792    if(DD[2*i-1]=="moverdim"){moverdim=DD[2*i];}
     5793  }
     5794  int nv=moverdim;
     5795  if(moverdim>nuv){moverdim=nuv;}
     5796
     5797  ideal vmov;
     5798  //string("T_nuv=",nuv,"; moverdim=",moverdim);
     5799  for(i=1;i<=moverdim;i++){vmov[size(vmov)+1]=var(i+nuv-moverdim);}
    60595800  if(size(GG)==0){GG=grobcov(F);}
    60605801  int j; int k; int te;
     
    60635804  list nGP;
    60645805  if (equalideals(B0,ideal(1)) )
    6065   {return(locus2(GG,F,DD));}
    6066     // if(version==2){return(locus2(GG,F,DD));}
    6067     //   else{return(locus0(GG,DD));}
    6068     // }
     5806  {return(locus2(GG,F,moverdim,vmov,na));}
    60695807  else
    60705808  {
    6071     n=nvars(RR);
    60725809    ideal vB;
    60735810    ideal N;
     
    60785815    attrib(N,"IsSB",1);
    60795816    N=std(N);
    6080    if((size(N))>=2)
     5817    if((size(N))>=2)
    60815818    {
    60825819       //def dN=dim(N);
     
    61205857      nGP=GG;
    61215858      " ";string("Unavoidable ",moverdim,"-dimensional locus");
    6122       //string("Try option 'vmov',ideal(of mover variables) to avoid some point of the mover");
    6123       //" ";"Elements of the basis of the generic segment in mover variables=";  N;" ";
    61245859      list L; return(L);
    61255860    }
    61265861  }
    6127   if(comment>0){"Input for locus2 GB="; nGP; "input for locus  F="; F;}
    6128   if(version==2)
    6129   {
    6130     def LL=locus2(nGP,F,DD);
    6131   }
    6132   else{def LL=locus0(nGP,DD);}
    6133   // kill @RP; kill @P; kill @R;
     5862
     5863//  if(comment>0){"Input for locus2 GB="; nGP; "input for locus  F="; F;}
     5864//  if(version==2)
     5865//  {
     5866//    "T_nGP enter for locus2="; nGP;
     5867//    def LL=locus2(nGP,F,moverdim,vmov,na);
     5868//  }
     5869//  else{ def LL=locus0(nGP,moverdim,vmov);  }
     5870
     5871  def LL=locus2(nGP,F,moverdim,vmov,na);
     5872
    61345873  return(LL);
    61355874}
    61365875example
    6137 {
    6138 "EXAMPLE:"; echo = 2;
     5876{ "EXAMPLE:"; echo = 2;
     5877
     5878// EXAMPLE 1
     5879
     5880// Conchoid, Pascal's Limacon.
     5881
     5882//  1. Given a circle: x1^2+y1^2-4
     5883//  2. and a mover point M(x1,y1) on it
     5884//  3. Consider the fix point P(0,2) on the circle
     5885//  4. Consider the line l passing through M and P
     5886//  5. The tracer T(x,y) are the points on l at fixed distance 1 to M.
     5887
    61395888if(defined(R)){kill R;}
    61405889ring R=(0,x,y),(x1,y1),dp;
     
    61465895locus(S96);
    61475896
    6148 // EXAMPLE
     5897// EXAMPLE 2
     5898
    61495899// Consider two parallel planes z1=-1 and z1=1, and two orthogonal parabolas on them.
    61505900// Determine the locus generated by the lines that rely the two parabolas
     
    61525902
    61535903if(defined(R)){kill R;}
    6154 ring R=(0,x,y,z),(lam,x2,y2,z2,x1,y1,z1), lp;
     5904ring R=(0,x,y,z),(x2,y2,z2,z1,y1,x1,lam), lp;
    61555905short=0;
    61565906
    61575907ideal L=z1+1,
    6158             x1^2-y1,
    6159             z2-1,
    6160             y2^2-x2,
    6161             4*x1*y2-1,
    6162             x-x1-lam*(x2-x1),
    6163             y-y1-lam*(y2-y1),
    6164             z-z1-lam*(z2-z1);
    6165 
    6166 locus(L);  // uses "moverdim",7
    6167 
    6168 // Now observe the effect on the antiimage of option moverdim:
    6169 // It does not take into account that lam is free.
    6170 locus(L,"moverdim",3);
    6171 
    6172 // Observe that with the option "mpverdim",3 the taxonomy becomes Special because considering only x1,y1,z1 as mover variables does not take into account that lam is free in the global anti-image.
     5908        x1^2-y1,
     5909        z2-1,
     5910        y2^2-x2,
     5911        4*x1*y2-1,
     5912        x-x1-lam*(x2-x1),
     5913        y-y1-lam*(y2-y1),
     5914        z-z1-lam*(z2-z1);
     5915
     5916locus(L);  // uses "moverdim",3
     5917// Observe the choose of the mover variables: the last 3 variables y1,x1,lam
     5918// If we choose x1,y1,z1 instead, the taxonomy becomes "Special" because
     5919// z1=-1 is fix and do not really correspond to the mover variables.
     5920
     5921// EXAMPLE 3 of parametric locus:
     5922
     5923// Determining the equation of a general ellipse;
     5924// Uncentered elipse;
     5925
     5926// Parameters  (a,b,a0,b0,p):
     5927//   a=large semiaxis, b=small semiaxis,
     5928//   (a0,b0) = center of the ellipse,
     5929//   (a1,b1) and (2*a0-a1,2*b0-b1) the focus,
     5930//   p the slope of the line of the a-axis of the ellipse.
     5931
     5932// Determine the equation of the ellipse.
     5933
     5934// We must use the option "numpar",5 in order to consider
     5935// the first 5 parameters as free parameters for the locus
     5936
     5937// Auxiliary variabes:
     5938//  d1=distance from focus (a1,b1) to the mover point M(x1,y1),
     5939//  d2=distance from focus (a2,b2) to the mover point M(x1,y1),
     5940//  f=focus distance= distance from (a0,b0) to (a1,b1).
     5941
     5942// Mover point (x1,y1) = tracer point (x,y).
     5943
     5944if(defined(R1)){kill R1;}
     5945ring R1=(0,a,b,a0,b0,p,x,y),(d1,d2,f,a1,b1,x1,y1),lp;
     5946
     5947ideal F3=b1-b0-p*(a1-a0),
     5948          //b2-b0+p*(a1-a0),
     5949          //a1+a2-2*a0,
     5950          //b1+b2-2*b0,
     5951          f^2-(a1-a0)^2-(b1-b0)^2,
     5952          f^2-a^2-b^2,
     5953          (x1-a1)^2+(y1-b1)^2-d1^2,
     5954          (x1-2*a0+a1)^2+(y1-2*b0+b1)^2-d2^2,
     5955          d1+d2-2*a,
     5956          x-x1,
     5957          y-y1;
     5958
     5959def G3=grobcov(F3);
     5960
     5961def Loc3=locus(F3,"grobcov",G3,"numpar",5); Loc3;
     5962
     5963// General ellipse:
     5964
     5965def C=Loc3[1][1][1];
     5966C;
     5967
     5968// Centered ellipse of semiaxes (a,b):
     5969
     5970def C0=subst(C,a0,0,b0,0,p,0);
     5971C0;
    61735972}
    61745973
     
    62116010}
    62126011example
    6213 {
    6214   "EXAMPLE:"; echo = 2;
     6012{ "EXAMPLE:"; echo = 2;
    62156013if(defined(R)){kill R;};
    62166014ring R=(0,a,b),(x,y),dp;
    62176015short=0;
     6016
    62186017// Concoid
    62196018ideal S96=x^2+y^2-4,(b-2)*x-a*y+2*a,(a-x)^2+(b-y)^2-1;
     
    63006099}
    63016100example
    6302 {
    6303   "EXAMPLE:"; echo = 2;
     6101{ "EXAMPLE:"; echo = 2;
    63046102if(defined(R)){kill R;}
    63056103ring R=(0,x,y),(x1,y1),dp;
    63066104short=0;
     6105
    63076106ideal S=x1^2+y1^2-4,(y-2)*x1-x*y1+2*x,(x-x1)^2+(y-y1)^2-1;
    63086107def L=locus(S);
     
    63296128       the parameteres of the hyper-surfaces, that are
    63306129       considered as variables of the parametric ring.
     6130       In the actual version, parametric envelope are accepted.
     6131       To include fixed parameters a1,..ap, to the problem, one must
     6132       declare them as the first parameters of the ring. if the
     6133       the number of free parameters is p, the option \"numpar\",p
     6134       is required.
    63316135       Calling sequence:
    6332        ring R=(0,x_1,..,x_n),(u_1,..,u_m),lp;
    6333        poly F=F(x_1,..,x_n,u_1,..,u_m);
    6334        ideal C=g_1(u_1,..u_m),..,g_s(u_1,..u_m);
     6136       ring R=(0,a1,..,ap,x_1,..,x_n),(u_1,..,u_m),lp;
     6137       poly F=F(a1,..ap,x_1,..,x_n,u_1,..,u_m);
     6138       ideal C=g_1(a1,..,ap,u_1,..u_m),..,g_s(a1,..ap,u_1,..u_m);
    63356139       envelop(F,C[,options]);   where s<m.
    6336        x1,..,xn are the trtacer variables.
     6140       x1,..,xn are the tracer variables.
    63376141       u_1,..,u_m are the auxiliary variables.
     6142       a1,..,ap are the fixed parameters if they exist
     6143       If the problem is a parametric envelope, and a's exist,
     6144       then the option \"numpar\",p  m must be given.
    63386145       By default the las n variables are the mover variables.
     6146       See the EXAMPLE of parametric envelop by calling
     6147       example envelop,
    63396148RETURN: The output is a list of the components [C_1, .. , C_n]
    63406149       of the locus. Each component is given by
     
    63656174       \"moverdim\", k: by default it is equal to n, the number of
    63666175       x-tracer variables.
     6176       \"numpar\",p  when fixed parameters are included
    63676177NOTE: grobcov and locus are called internally.
    6368        The basering R, must be of the form Q[x][u]
    6369        (x=parameters, u=variables).
     6178       The basering R, must be of the form Q[a,x][u]
     6179       (x=variables, u=auxiliary variables), (a fixed parameters).
    63706180       This routine uses the generalized definition of envelop
    63716181       introduced in the book
     
    63766186{
    63776187  def RR=basering;
    6378   int tes=0; int i;   int j;  int k; int m;
    6379   int d;
    6380   int dp;
    6381   ideal BBB;
    6382   // new
    6383   ideal CC=C;
    6384   CC[size(CC)+1]=F;
    6385   int movreal=size(variables(CC));
    6386   //Options
     6188  list LRR=ringlist(RR);
     6189  int nax=size(LRR[1][2]);
     6190  int nuv=size(LRR[2]);
     6191
    63876192  list DD=#;
    6388 
    6389   list Rx=ringlist(RR);
    6390   ideal vmov;
    6391   int moverdim=movreal;         //size(Rx[1][2]);
    6392   int nv=nvars(RR);
    6393   if(moverdim>nv){moverdim=nv;}
    6394   for(i=1;i<=moverdim;i++){vmov[size(vmov)+1]=var(i+nv-moverdim);}
    6395   //for(i=1;i<=nv;i++)  {vmov[size(vmov)+1]=var(i);}
    6396   int version=2;
    6397   int comment=0;
    6398   ideal Fm;
    6399   int tax=1;
    6400 
    6401   int familyinfo;
    6402   for(i=1;i<=(size(DD) div 2);i++)
    6403   {
    6404     if(DD[2*i-1]=="anti-image"){tax=DD[2*i];}
    6405     if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    6406     if(DD[2*i-1]=="version"){version=2;DD[2*i]=2;}
    6407     if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    6408     if(DD[2*i-1]=="familyinfo"){familyinfo=0;DD[2*i]=0;}
    6409   };
    6410   int ng=size(C);
    6411   ideal S=F;
    6412   for(i=1;i<=size(C);i++){S[size(S)+1]=C[i];}
    6413   int s=nv-ng;
    6414   if(s>0)
    6415   {
    6416     matrix M[ng+1][ng+1];
    6417     def cc=comb(nv,ng+1);
     6193  int na=0;
     6194  int nx=nax;
     6195  int nu=0;
     6196  int nv=nuv;
     6197  int i; int j; int k;
     6198  //string("T_ nax=",nax,"; nx=",nx,"; nuv=",nuv,"; nv=",nv);
     6199  int tnumpar=0;
     6200  // int tnumvar=0;
     6201  //"T_DD="; DD;
     6202  for(i=1;i<=size(DD) div 2;i++)
     6203  {
     6204    if(DD[2*i-1]=="numpar"){na=DD[2*i];tnumpar=1;}
     6205    // if(DD[2*i-1]=="numvar"){nv=DD[2*i];tnumvar=1;}
     6206  }
     6207  if(tnumpar==0){DD[size(DD)+1]="numpar";  DD[size(DD)+1]=na;}
     6208  // if(tnumvar==0){DD[size(DD)+1]="numvar";  DD[size(DD)+1]=nv;}
     6209  nx=nax-na;
     6210  nu=nuv-nv;
     6211  //string("T_ nax=",nax,"; nx=",nx,"; nuv=",nuv,"; nv=",nv);
     6212  ideal Vnv;
     6213  ideal Vnonv;
     6214  for(i=1;i<=nu;i++){Vnonv[size(Vnonv)+1]=var(i);}
     6215  //"T_Vnonv="; Vnonv;
     6216  for(i=nu+1;i<=nuv;i++){Vnv[size(Vnv)+1]=var(i);}
     6217  //"T_Vnv="; Vnv;
     6218  ideal Cnor;
     6219  ideal Cr=F;
     6220  for(i=1;i<=size(C);i++)
     6221  {
     6222    if(subset(variables(C[i]),Vnonv)){Cnor[size(Cnor)+1]=C[i];}
     6223    else{Cr[size(Cr)+1]=C[i];}
     6224  }
     6225  int nr=size(Cr);
     6226  //string("T_nr=", nr,"; nv=",nv);
     6227  if(nr>0)
     6228  {
     6229    matrix M[nr][nr];
     6230    def cc=comb(nv,nr);
     6231    //"T_cc="; cc;
     6232    //string("T_nv=",nv," nr=",nr);
    64186233    poly J;
    64196234    for(k=1;k<=size(cc);k++)
    64206235    {
    6421       for(j=1;j<=ng+1;j++)
     6236      for(i=1;i<=nr;i++)
    64226237      {
    6423         M[1,j]=diff(F,var(cc[k][j]));
    6424       }
    6425       for(i=1;i<=ng;i++)
    6426       {
    6427         for(j=1;j<=ng+1;j++)
     6238        for(j=1;j<=nr;j++)
    64286239        {
    6429           M[i+1,j]=diff(C[i],var(cc[k][j]));
     6240          M[i,j]=diff(Cr[i],var(cc[k][j]));
    64306241        }
    64316242      }
    64326243      J=det(M);
    6433       S[size(S)+1]=J;
    6434     }
    6435  }
    6436  if(comment>0){"System S before grobcov ="; S;}
    6437   //def G=grobcov(S,DD);
    6438   list HHH;
    6439   DD[size(DD)+1]="moverdim";   // new option
    6440   DD[size(DD)+1]=moverdim;
    6441   // string("moverdim=",moverdim);
     6244      Cr[size(Cr)+1]=J;
     6245    }
     6246  }
     6247  ideal S=Cnor;
     6248  for(i=1;i<=size(C);i++){S[size(S)+1]=C[i];}
     6249  for(i=1;i<=size(Cr);i++){S[size(S)+1]=Cr[i];}
     6250  //"T_S="; S;
    64426251  def L=locus(S,DD);
    6443   // "L="; L;
    6444   list GL;
    6445   ideal fam; ideal env;
    6446 
    6447   def P=ring(Rx[1]);
    6448   Rx[1]=0;
    6449   def D=ring(Rx);
    6450   def RP=P+D;
    6451   list LL;
    6452   list NormalComp;
    6453   ideal Gi;
    6454   ideal BBBB;
    6455   poly B0;
    64566252  return(L);
    64576253}
    64586254example
    6459 {
    6460   "EXAMPLE:"; echo = 2;
    6461 if(defined(R)){kill R;}
    6462 ring R=(0,x,y),(r,s,y1,x1),lp;
    6463 poly F=(x-x1)^2+(y-y1)^2-r;
    6464 ideal g=(x1-2*(s+r))^2+(y1-s)^2-s;
    6465 
    6466 def E=envelop(F,g);
    6467 E;
    6468 
    6469 E[1][1];
    6470 
    6471 string("Taxonomy=", E[1][3][2]);
    6472 
    6473 
    6474 // EXAMPLE
     6255{ "EAXMPLE:"; echo=2;
     6256
     6257// EXAMPLE 1
    64756258// Steiner Deltoid
    64766259// 1. Consider the circle x1^2+y1^2-1=0, and a mover point M(x1,y1) on it.
    64776260// 2. Consider the triangle A(0,1), B(-1,0), C(1,0).
    64786261// 3. Consider lines passing through M perpendicular to two sides of ABC triangle.
    6479 // 4. Obtain the envelop of the lines above.
     6262// 4. Determine the envelope of the lines above.
     6263
    64806264if(defined(R)){kill R;}
    64816265ring R=(0,x,y),(x1,y1,x2,y2),lp;
    64826266short=0;
     6267
    64836268ideal C=(x1)^2+(y1)^2-1,
    64846269             x2+y2-1,
     
    64866271matrix M[3][3]=x,y,1,x2,y2,1,x1,0,1;
    64876272poly F=det(M);
    6488 // Curves Family F
     6273
     6274// The lines of family F are
    64896275F;
    6490 // Conditions C=
     6276
     6277// The conditions C are
    64916278C;
     6279
    64926280envelop(F,C);
     6281
     6282// EXAMPLE 2
     6283// Parametric envelope
     6284
     6285// Let c be the circle centered at the origin O(0,0) and having radius 1.
     6286// M(x1,y1) be a mover point gliding on c.
     6287// Let A(a0,b0) be a parametric fixed point:
     6288// Consider the set of lines parallel to the line AO passing thoug M.
     6289
     6290// Determine the envelope of these lines
     6291
     6292// We let the fixed point A coordinates as free parameters of the envelope.
     6293// We have to declare the existence of two parameters when
     6294// defining the ring in which we call envelop,
     6295// and set a0,b0 as the first variables of the parametric ring
     6296// The ring is thus
     6297
     6298if(defined(R1)){kill R1;}
     6299ring R1=(0,a0,b0,x,y),(x1,y1),lp;
     6300short=0;
     6301
     6302// The lines are  F1
     6303poly F1=b0*(x-x1)-a0*(y-y1);
     6304
     6305// and the mover is on the circle c
     6306ideal C1=x1^2+y1^2-1;
     6307// The call is thus
     6308
     6309def E1=envelop(F1,C1,"numpar",2);
     6310E1;
     6311
     6312// The interesting first component  EC1 is
     6313def EC1=E1[1][1][1];
     6314EC1;
     6315
     6316// that is equivalent to  (a0*y-b0*x)^2-a0^2-b0^2.
     6317// As expected it consists of the two lines
     6318//    a0*y-b0*x - sqrt(a0^2+b0^2),
     6319//    a0*y-b0*x + sqrt(a0^2+b0^2),
     6320// parallel to the line OM passing at the
     6321// points of the circle in the line perpendicular to OA.
     6322
     6323// EXAMPLE 3
     6324// Parametric envelope
     6325
     6326// Let c be the circle centered at the origin O(a1,b1) and having radiusr,
     6327// where a1,b1,r are fixed parameters
     6328// M(x1,y1) be a mover point gliding on c.
     6329// Let A(a0,b0) be a parametric fixed point:
     6330// Consider the set of lines parallel to the line AO passing thoug M.
     6331
     6332// Determine the envelope of these lines
     6333
     6334// We let the fixed point A,point M and r as free parameters of the envelope.
     6335// We have to declare the existence of 5 parameters when
     6336// defining the ring in which we call envelop,
     6337// and set a0,b0,a1,b1,r as the first variables of the parametric ring
     6338// The ring is thus
     6339
     6340if(defined(R1)){kill R1;}
     6341ring R1=(0,a0,b0,a1,b1,r,x,y),(x1,y1),lp;
     6342short=0;
     6343
     6344// The lines are  F1
     6345poly F1=b0*(x-x1)-a0*(y-y1);
     6346
     6347// and the mover is on the circle c
     6348ideal C1=(x1-a1)^2+(y1-b1)^2-r^2;
     6349// The call is thus
     6350
     6351def E1=envelop(F1,C1,"numpar",5);
     6352E1;
     6353
     6354// The interesting first component  EC1 is
     6355def EC1=E1[1][1][1];
     6356EC1;
     6357
     6358// which corresponds to the product of two lines
     6359// parallel to the line AM and intercepting the circle
     6360// on the intersection of the line perpendicuar
     6361// to line AM passing through A
    64936362}
    64946363
     
    65386407  ideal BBB;
    65396408  //Options
    6540   list DD=#;
     6409//    list DD=#;
    65416410  ideal vmov;
    65426411  int nv=nvars(RR);
    65436412  for(i=1;i<=nv;i++){vmov[size(vmov)+1]=var(i);}
    6544   int numpars=size(ringlist(RR)[1][2]);
    6545   int version=0;
    6546   if(nv<4){version=1;}
     6413//  int numpars=size(ringlist(RR)[1][2]);
     6414//  int version=0;
     6415//  if(nv<4){version=1;}
    65476416  int comment=0;
    6548   int familyinfo;
     6417  int familyinfo=0;
    65496418  ideal Fm;
    6550   for(i=1;i<=(size(DD) div 2);i++)
    6551   {
    6552     if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    6553 //    if(DD[2*i-1]=="version"){version=DD[2*i];}
    6554     if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    6555     if(DD[2*i-1]=="familyinfo"){familyinfo=DD[2*i];}
    6556     if(DD[2*i-1]=="moreinfo"){moreinfo=DD[2*i];}
    6557   };
    6558   DD=list("vmov",vmov,"version",version,"comment",comment);
     6419//    for(i=1;i<=(size(DD) div 2);i++)
     6420//    {
     6421//      if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
     6422//  //    if(DD[2*i-1]=="version"){version=DD[2*i];}
     6423//      if(DD[2*i-1]=="comment"){comment=DD[2*i];}
     6424//      if(DD[2*i-1]=="familyinfo"){familyinfo=DD[2*i];}
     6425//      if(DD[2*i-1]=="moreinfo"){moreinfo=DD[2*i];}
     6426//    };
     6427//    DD=list("vmov",vmov,"comment",comment); // ,"version",version
    65596428  int ng=size(C);
    65606429  ideal S=F;
     
    65886457 }
    65896458 //if(comment>0){"System S before grobcov ="; S;}
    6590   def G=grobcov(S,DD);
     6459 //"T_S="; S;
     6460  def G=grobcov(S); // ,DD
     6461  //"T_G=";G;
    65916462  list GG;
    65926463  for(i=2;i<=size(G);i++)
     
    65956466  }
    65966467  G=GG;
     6468  //"T_G=";G;
    65976469  if(moreinfo>0){return(G);}
    65986470  else
     
    66056477      //string("T_G[",i,"][3][1][1][1]="); G[i][3][1][1][1];
    66066478      //string("T_EE="); EE;
    6607       if(G[i][3][1][1][1]==EE)
     6479      if(equalideals(G[i][3][1][1],EE))
    66086480      {
    66096481         t=1;
     
    66176489}
    66186490example
    6619 {
    6620   "EXAMPLE:"; echo = 2;
    6621 
     6491{ "EXAMPLE:"; echo = 2;
    66226492if(defined(R)){kill R;}
    66236493ring R=(0,x,y),(r,s,y1,x1),lp;
     
    66476517
    66486518proc FamElemsAtEnvCompPoints(poly F,ideal C, ideal E,list #)
    6649 "USAGE: FamElemsAtEnvCompPoints(poly F,ideal C,ideal E);
     6519"USAGE: FamElemsAtEnvCompPoints(poly F,ideal C,poly E);
    66506520       poly F must be the family of hyper-surfaces whose
    66516521       envelope is analyzed. It must be defined in the ring
     
    66916561  int i;
    66926562  int moreinfo=0;
     6563  int familyinfo=0;
     6564  int comment=0;
     6565  int numpar=0;
     6566  ideal vmov;
    66936567  list DD=#;
    66946568  for(i=1;i<=(size(DD) div 2);i++)
    66956569  {
    66966570    if(DD[2*i-1]=="vmov"){vmov=DD[2*i];}
    6697 //    if(DD[2*i-1]=="version"){version=DD[2*i];}
    66986571    if(DD[2*i-1]=="comment"){comment=DD[2*i];}
    66996572    if(DD[2*i-1]=="familyinfo"){familyinfo=DD[2*i];}
    67006573    if(DD[2*i-1]=="moreinfo"){moreinfo=DD[2*i];}
     6574    if(DD[2*i-1]=="numpar"){numpar=DD[2*i];}
    67016575  };
    67026576  ideal S=C;
     
    67236597      //string("T_G[",i,"][3][1][1][1]="); G[i][3][1][1][1];
    67246598      //string("T_EE="); EE;
    6725       if(G[i][3][1][1][1]==EE)
     6599      if(G[i][3][1][1][1]==E)
    67266600      {
    67276601         t=1;
     
    67346608}
    67356609example
    6736 {
    6737   "EXAMPLE:"; echo = 2;
    6738 
    6739 if(defined(R)){kill R;}
    6740 ring R=(0,x,y),(t),dp;
    6741 short=0;
    6742 poly F=(x-5*t)^2+y^2-9*t^2;
    6743 ideal C;
    6744 
    6745 def Env=envelop(F,C);
    6746 Env;
     6610{ "EXAMPLE:"; echo = 2;
     6611 if(defined(R)){kill R;}
     6612 ring R=(0,x,y),(t),dp;
     6613 short=0;
     6614 poly F=(x-5*t)^2+y^2-9*t^2;
     6615 ideal C;
     6616
     6617 def Env=envelop(F,C);
     6618 Env;
    67476619
    67486620// E is a component of the envelope:
    6749 ideal E=Env[1][1];
    6750 
    6751 E;
    6752 
    6753 def A=AssocTanToEnv(F,C,E);
    6754 A;
     6621 def E=Env[1][1][1];
     6622 E;
     6623
     6624 def A=AssocTanToEnv(F,C,E);
     6625 A;
    67556626
    67566627// The basis of the parameter values of the associated
     
    67616632//    element at (x,y) is
    67626633
    6763 subst(F,t,-(5/12)*y);
    6764 
    6765 def FE=FamElemsAtEnvCompPoints(F,C,E);
    6766 FE;
    6767 
    6768 factorize(FE[2][1]);
     6634 subst(F,t,-(5/12)*y);
     6635
     6636 def FE=FamElemsAtEnvCompPoints(F,C,E);
     6637 FE;
     6638
     6639 factorize(FE[2][1]);
    67696640
    67706641// Thus the unique family element passing through the envelope point (x,y)
     
    67726643
    67736644// EXAMPLE:
    6774 if(defined(R)){kill R;}
    6775 ring R=(0,x,y),(r,s,y1,x1),lp;
    6776 
    6777 poly F=(x-x1)^2+(y-y1)^2-r;
    6778 ideal g=(x1-2*(s+r))^2+(y1-s)^2-s;
    6779 
    6780 def E=envelop(F,g);
    6781 E;
    6782 
    6783 def A=AssocTanToEnv(F,g,E[1][1][1]);
    6784 A;
    6785 
    6786 def M1=coef(A[2][1],x1);
    6787 def M2=coef(A[2][2],y1);
    6788 def M3=coef(A[2][3],s);
    6789 def M4=coef(A[2][4],r);
     6645 if(defined(R)){kill R;}
     6646 ring R=(0,x,y),(r,s,y1,x1),lp;
     6647
     6648 poly F=(x-x1)^2+(y-y1)^2-r;
     6649 ideal g=(x1-2*(s+r))^2+(y1-s)^2-s;
     6650
     6651 def E=envelop(F,g);
     6652 E;
     6653
     6654 def A=AssocTanToEnv(F,g,E[1][1][1]);
     6655 A;
     6656
     6657 def M1=coef(A[2][1],x1);
     6658 def M2=coef(A[2][2],y1);
     6659 def M3=coef(A[2][3],s);
     6660 def M4=coef(A[2][4],r);
    67906661
    67916662// The parameter values corresponding to the family
    67926663//    element tangent at point (x,y) of the envelope are:
    6793 "x1=";-M1[2,2]/M1[2,1];
    6794 
    6795 "y1=";-M2[2,2]/M2[2,1];
    6796 
    6797 "s=";-M3[2,2]/M3[2,1];
    6798 
    6799 "r=";-M4[2,2]/M4[2,1];
     6664 "x1=";-M1[2,2]/M1[2,1];
     6665
     6666 "y1=";-M2[2,2]/M2[2,1];
     6667
     6668 "s=";-M3[2,2]/M3[2,1];
     6669
     6670 "r=";-M4[2,2]/M4[2,1];
    68006671
    68016672// Now detect if there are other family elements passing at this point:
    6802 def FE=FamElemsAtEnvCompPoints(F,g,E[1][1][1]);
    6803 FE;
     6673 def FE=FamElemsAtEnvCompPoints(F,g,E[1][1][1]);
     6674 FE;
    68046675
    68056676// FE[1] is the set of lpp. It has dimension 4-2=2.
     
    68526723}
    68536724example
    6854 {
    6855   "EXAMPLE:"; echo = 2;
     6725{ "EXAMPLE:"; echo = 2;
    68566726if(defined(R)){kill R;}
    68576727ring R=(0,a,b,c),(x,y),dp;
     
    68626732}
    68636733
    6864 // AddLocus: auxiliary routine for locus0 that computes the components of the constructible:
     6734// AddLocus: auxilliary routine for locus0 that computes the components of the constructible:
    68656735// Input:  the list of locally closed sets to be added, each with its type as third argument
    68666736//     L=[ [LC[11],..,LC[1k_1],.., [LC[r1],..,LC[rk_r] ] where
     
    71286998    }
    71296999     // B=reduced basis on CR
     7000     //"T_PR="; PR;
     7001     //"T_CR="; CR;
     7002     //"T_B="; B;
    71307003    if(rep==1){return(list(lcc,B,CR));}
    71317004    else{return(list(lcc,B,PR));}
     
    71427015}
    71437016example
    7144 {
    7145   "EXAMPLE:"; echo = 2;
    7146   if(defined(RE)){kill RE;}
    7147   ring RE=(0,a,b,c,d,e,f),(x,y),lp;
    7148   ideal F=a*x^2+b*x*y+c*y^2,d*x^2+e*x*y+f*y^2;
    7149   ideal A=a*e-b*d;
     7017{ "EXAMPLE:"; echo = 2;
     7018if(defined(RE)){kill RE;}
     7019ring RE=(0,a,b,c,d,e,f),(x,y),lp;
     7020ideal F=a*x^2+b*x*y+c*y^2,d*x^2+e*x*y+f*y^2;
     7021ideal A=a*e-b*d;
    71507022
    71517023WLemma(F,A);
     
    72287100    {
    72297101      G[size(G)+1]=G0;
     7102      //"T_G0="; G0;
    72307103      N=G0[3][1][2];
     7104      //"T_N="; N;
    72317105      for(i=1;i<=size(N);i++)
    72327106      {
     
    72367110          Epend[size(Epend)+1]=N[i];
    72377111        }
     7112        //"T_i="; i;
    72387113      }
     7114      //"T_Etot="; Etot;
     7115      //"T_Epend=";Epend;
    72397116    }
    72407117  }
     
    72427119}
    72437120example
    7244 {
    7245   "EXAMPLE:"; echo = 2;
    7246   if(defined(RRR)){kill RRR;}
    7247   ring RRR=(0,b,c,d,e,f),(x,y,t),lp;
    7248   short=0;
    7249   ideal S=x^2+2*c*x*y+2*d*x*t+b*y^2+2*e*y*t+f*t^2,
     7121{ "EXAMPLE:"; echo = 2;
     7122 if(defined(RRR)){kill RRR;}
     7123 ring RRR=(0,b,c,d,e,f),(x,y,t),lp;
     7124 short=0;
     7125 ideal S=x^2+2*c*x*y+2*d*x*t+b*y^2+2*e*y*t+f*t^2,
    72507126            x+c*y+d*t,c*x+b*y+e*t;
    7251   grobcov(S);
    7252   WLcgs(S);
     7127 grobcov(S);
     7128 WLcgs(S);
    72537129}
    72547130
     
    73127188  int i; int j;
    73137189  ideal J=idminusid(E,N);
     7190  //"T_J="; J;
    73147191  for(i=1;i<=size(F);i++)
    73157192  {
     
    74227299}
    74237300example
    7424 {
    7425 "EXAMPLE:"; echo = 2;
    7426   if (defined(R)) {kill R;}
    7427   ring R=(0,x,y),(x1,y1,x2,y2),lp;
    7428   ideal F=-y*x1+(x-1)*y1+y,
    7429             (x-1)*(x1+1)+y*y1,
    7430             -y*x2+(x+1)*y2-y,
    7431             (x+1)*(x2-1)+y*y2,
    7432             (x1-x)^2+y1^2-(x1-x)^2-y2^2;
     7301{ "EAXMPLE:"; echo = 2;
     7302if (defined(R)) {kill R;}
     7303ring R=(0,x,y),(x1,y1,x2,y2),lp;
     7304ideal F=-y*x1+(x-1)*y1+y,
     7305         (x-1)*(x1+1)+y*y1,
     7306         -y*x2+(x+1)*y2-y,
     7307         (x+1)*(x2-1)+y*y2,
     7308         (x1-x)^2+y1^2-(x1-x)^2-y2^2;
    74337309
    74347310def G=grobcov(F,"rep",1);
    7435 
    74367311G;
    74377312
    74387313def L=Grob1Levels(G);
    7439 
    74407314L;
    74417315
    74427316def LL=Levels(L);
    7443 
    74447317LL;
    74457318}
     
    74637336  setring(RP);
    74647337  def SP=imap(RR,S);
     7338  //"T_SP="; SP;
     7339  //"T_typeof(SP[1])="; typeof(SP[1]);
    74657340  ideal EP;
    74667341  EP=SP[1];
     
    74767351}
    74777352example
    7478 {
    7479 "EXAMPLE:"; echo = 2;
     7353{ "EAXMPLE:"; echo = 2;
    74807354if(defined(R)){kill R;}
    74817355ring R=(0,x,y,z),(x1,y1,z1),lp;
     7356
    74827357ideal I1=x+y*z*x1;
    74837358ideal I2=x-y*z*y1;
     
    74887363intersectpar(S);
    74897364}
    7490 
    74917365
    74927366proc ADGT(ideal H,ideal T,ideal H1,ideal T1,list #)
     
    76147488    int r=m div 2;
    76157489    if(m mod 2==0){r=r-1;}
    7616     for(i=1;i<=r;i++)
     7490    //"L0="; L0;
     7491     for(i=1;i<=r;i++)
    76177492    {
    76187493      G1=grobcov(F,"null",L0[2*i],"nonnull",L0[2*i+1],"rep",1);
     
    76777552}
    76787553example
    7679 {
    7680 "EXAMPLE:"; echo = 2;
    7681 
     7554{ "EXAMPLE:"; echo = 2;
    76827555// Determine the supplementary conditions
    76837556// for the non-degenerate  triangle A(-1,0), B(1,0), C(x,y)
     
    76917564
    76927565ideal H=-y*x1+(x-1)*y1+y,
    7693          (x-1)*(x1+1)+y*y1,
    7694          -y*x2+(x+1)*y2-y,
    7695          (x+1)*(x2-1)+y*y2;
     7566        (x-1)*(x1+1)+y*y1,
     7567        -y*x2+(x+1)*y2-y,
     7568        (x+1)*(x2-1)+y*y2;
    76967569
    76977570// Thesis T: the orthic triangle is isosceles
    7698 ideal T=(x1-x)^2+y1^2-(x2-x)^2-y2^2;
     7571 ideal T=(x1-x)^2+y1^2-(x2-x)^2-y2^2;
    76997572
    77007573// Negative Hypothesis H1: ABC is non-degenerate
    77017574ideal H1=y;
    77027575
    7703  // Negative Thesis T1: the orthic triangle is non-degenerate
     7576// Negative Thesis T1: the orthic triangle is non-degenerate
    77047577ideal T1=x*(y1-y2)-y*(x1-x2)+x1*y2-x2*y1;
    77057578
     
    77107583ADGT(H,T,H1,T1);
    77117584
    7712 // Now using difference of constructible sets for negative hypothesis and thesis
     7585// Now using diference of constructible sets for negative hypothesis and thesis
    77137586
    77147587ADGT(H,T,H1,T1,"neg",1);
    77157588
    77167589// The results are identical using both methods for the negative propositions
    7717 //  - Rabinovitch or
    7718 //  - DifConsLCSets
    7719 
    7720 // EXAMPLE
     7590// - Rabinovitch or
     7591// - DifConsLCSets
     7592
     7593// EXAMPLE 2
    77217594
    77227595// Automatic Theorem Proving.
     
    77417614
    77427615if (defined(R1)){kill R1;}
    7743 
    77447616ring R1=(0,a,b),(x1,y1,x2,y2,x0,y0,x3,y3,r2),dp;
    7745 
    77467617short=0;
    77477618
     
    77637634       (x0+a-2*x3)^2+(y0+b-2*y3)^2-r2;
    77647635
    7765  ADGT(H,T,b,1);
     7636ADGT(H,T,b,1);
    77667637
    77677638// Thus the nine point circle theorem is true for all real points excepting b=0.
     
    77977668    list G2=grobcov(F,"rep",1);
    77987669    def L2=Grob1Levels(G2);
     7670    //"L2="; L2;
    77997671    ideal FN=T1;
    78007672    if(not(equalideals(H1,ideal(1))))
     
    78117683      list G3=grobcov(F,"rep",1);
    78127684      def L3=Grob1Levels(G3);
     7685      //"L3="; L3;
    78137686    }
    78147687    def LL=DifConsLCSets(L2,L3);
     7688    //"LL="; LL;
    78157689    LL=ConsLevels(LL);
    78167690    LL=Levels(LL);
     7691     //"LL="; LL;
    78177692    if(rep==1){return(LL);}
    78187693    else
Note: See TracChangeset for help on using the changeset viewer.