Changeset b1fb09 in git for Singular/LIB/finvar.lib


Ignore:
Timestamp:
Jan 23, 2007, 2:03:59 PM (17 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', '17f1d200f27c5bd38f5dfc6e8a0879242279d1d8')
Children:
a2c203153a35ac48102771cb18fc358aab446969
Parents:
a38082faca14f3369605c3c541368da3a8a95d58
Message:
*hannes: without secondary_charp_degbound


git-svn-id: file:///usr/local/Singular/svn/trunk@9745 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/finvar.lib

    ra38082 rb1fb09  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: finvar.lib,v 1.53 2007-01-22 16:34:33 Singular Exp $"
     2version="$Id: finvar.lib,v 1.54 2007-01-23 13:03:59 Singular Exp $"
    33category="Invariant theory";
    44info="
    … …  
    4040 irred_secondary_char0()           irreducible secondary invariants in char 0
    4141 secondary_charp()                 secondary invariants in char p
    42  secondary_charp_degbound()        finds all (irred.) secondary invariants in char p
    43                                    up to a given degree bound, without using Molien
    44                                    series or Reynolds operator
    4542 secondary_no_molien()             secondary invariants, without Molien series
    4643                                   but with Reynolds operator
    … …  
    16211618           print(invariant_basis(2,A));
    16221619}
    1623 ///////////////////////////////////////////////////////////////////////////////
    1624 
    1625 static proc invariant_basis_modS (int g,ideal sS, list #)
    1626 "USAGE:   invariant_basis(g,sP,G1,G2,...);
    1627          g: an <int> indicating of which degree (>0) the homogeneous basis
    1628          shoud be,
    1629          sS: normal forms of subsequently found sec. inv. in degree g,
    1630          G1,G2,...: <matrices> generating a finite matrix group
    1631 RETURNS: invariants of degree g
    1632 THEORY:  A general polynomial of degree g is generated and the generators of
    1633          the matrix group applied. The difference ought to be 0 and this way a
    1634          system of linear equations is created. It is solved by computing
    1635          syzygies.
    1636 "
    1637 { if (g<=0)
    1638   { "ERROR:   the first parameter should be > 0";
    1639     return();
    1640   }
    1641   def br=basering;
    1642   ideal mon=sort(kbase(sS,g))[1];      // needed for constructing a general
    1643   int m=ncols(mon);                    // homogeneous polynomial of degree g
    1644   mon=sort(mon,intvec(m..1))[1];
    1645   int a=size(#);
    1646   int i;
    1647   int n=nvars(br);
    1648  //---------------------- checking that the input is ok -----------------------
    1649   for (i=1;i<=a;i++)
    1650   { if (typeof(#[i])=="matrix")
    1651     { if (nrows(#[i])==n && ncols(#[i])==n)
    1652       { matrix G(i)=#[i];
    1653       }
    1654       else
    1655       { "ERROR:   the number of variables of the base ring needs to be the same";
    1656         "         as the dimension of the square matrices";
    1657         return();
    1658       }
    1659     }
    1660     else
    1661     { "ERROR:   the last parameters should be a list of matrices";
    1662       return();
    1663     }
    1664   }
    1665  //----------------------------------------------------------------------------
    1666   execute("ring T=("+charstr(br)+"),("+varstr(br)+",p(1..m)),lp;");
    1667   // p(1..m) are the general coefficients of the general polynomial of degree g
    1668   execute("ideal vars="+varstr(br)+";");
    1669   map f;
    1670   ideal mon=imap(br,mon);
    1671   poly P=0;
    1672   for (i=m;i>=1;i--)
    1673   { P=P+p(i)*mon[i];                   // P is the general polynomial
    1674   }
    1675   ideal I;                             // will help substituting variables in P
    1676                                        // by linear combinations of variables -
    1677   poly Pnew,temp;                      // Pnew is P with substitutions -
    1678   matrix S[m*a][m];                    // will contain system of linear
    1679                                        // equations
    1680   int j,k;
    1681  //------------------- building the system of linear equations ----------------
    1682   for (i=1;i<=a;i++)
    1683   { I=ideal(matrix(vars)*transpose(imap(br,G(i))));
    1684     I=I,p(1..m);
    1685     f=T,I;
    1686     Pnew=f(P);
    1687     for (j=1;j<=m;j++)
    1688     { temp=P/mon[j]-Pnew/mon[j];
    1689       for (k=1;k<=m;k++)
    1690       { S[m*(i-1)+j,k]=temp/p(k);
    1691       }
    1692     }
    1693   }
    1694  //----------------------------------------------------------------------------
    1695   setring br;
    1696   map f=T,ideal(0);
    1697   matrix S=f(S);
    1698   matrix s=matrix(syz(S));             // s contains a basis of the space of
    1699                                        // solutions -
    1700   ideal I=ideal(matrix(mon)*s);        // I contains a basis of homogeneous
    1701   if (I[1]<>0)                         // invariants of degree d
    1702   { for (i=1;i<=ncols(I);i++)
    1703     { I[i]=I[i]/leadcoef(I[i]);        // setting leading coefficients to 1
    1704     }
    1705   }
    1706   return(I);
    1707 }
    1708 
    17091620///////////////////////////////////////////////////////////////////////////////
    17101621
    … …  
    56415552 //---------------- checking input and setting verbose mode -------------------
    56425553  if (char(br)==0)
    5643   { "ERROR:   secondary_charp should only be used with rings of characteristic p>0.";
    5644     return();
     5554  { ERROR("ERROR:   secondary_charp should only be used with rings of characteristic p>0.");
    56455555  }
    56465556  int i;
    … …  
    56695579      }
    56705580      else
    5671       { "ERROR:   If the last optional parameter is a string, it should be \"old\".";
    5672         return(matrix(ideal()),matrix(ideal()));
     5581      { ERROR("ERROR:   If the last optional parameter is a string, it should be \"old\".");
     5582        //return(matrix(ideal()),matrix(ideal()));
    56735583      }
    56745584    }
    … …  
    56795589                                       // we should get
    56805590  if (ncols(P)<>n)
    5681   { "ERROR:   The first parameter ought to be the matrix of the primary";
    5682     "         invariants."
    5683     return();
     5591  { ERROR("ERROR:   The first parameter ought to be the matrix of the primary"+newline+
     5592    "         invariants.");
    56845593  }
    56855594  if (ncols(REY)<>n)
    5686   { "ERROR:   The second parameter ought to be the Reynolds operator."
    5687     return();
     5595  { ERROR("ERROR:   The second parameter ought to be the Reynolds operator.");
    56885596  }
    56895597  if (typeof(`ring_name`)<>"ring")
    5690   { "ERROR:   The <string> should give the name of the ring where the Molien."
    5691     "         series is stored.";
    5692     return();
     5598  { ERROR("ERROR:   The <string> should give the name of the ring where the Molien."+newline+
     5599    "         series is stored.");
    56935600  }
    56945601  if (v && voice==2)
    … …  
    60825989        }
    60835990        else
    6084         { "ERROR:   the third parameter should be an <intvec>";
    6085           return();
     5991        { ERROR("ERROR:   the third parameter should be an <intvec>");
    60865992        }
    60875993      }
    … …  
    60956001        }
    60966002        else
    6097         { "ERROR:   the third parameter should be an <intvec>";
    6098           return();
     6003        { ERROR("ERROR:   the third parameter should be an <intvec>");
    60996004        }
    61006005      }
    61016006      else
    6102       { "ERROR:   wrong list of parameters";
    6103         return();
     6007      { ERROR("ERROR:   wrong list of parameters");
    61046008      }
    61056009    }
    … …  
    61076011  else
    61086012  { if (size(#)>2)
    6109     { "ERROR:   there are too many parameters";
    6110       return();
     6013    { ERROR("ERROR:   there are too many parameters");
    61116014    }
    61126015    int v=0;
    … …  
    61176020                                       // we should get
    61186021  if (ncols(P)<>n)
    6119   { "ERROR:   The first parameter ought to be the matrix of the primary";
    6120     "         invariants."
    6121     return();
     6022  { ERROR("ERROR:   The first parameter ought to be the matrix of the primary"+newline+
     6023    "         invariants.");
    61226024  }
    61236025  if (ncols(REY)<>n)
    6124   { "ERROR:   The second parameter ought to be the Reynolds operator."
    6125     return();
     6026  { ERROR("ERROR:   The second parameter ought to be the Reynolds operator.");
    61266027  }
    61276028  if (v && voice==2)
    … …  
    62356136        }
    62366137        else
    6237         { "ERROR:   the third parameter should be an <intvec>";
    6238           return();
     6138        { ERROR("ERROR:   the third parameter should be an <intvec>");
    62396139        }
    62406140      }
    … …  
    62486148        }
    62496149        else
    6250         { "ERROR:   the third parameter should be an <intvec>";
    6251           return();
     6150        { ERROR("ERROR:   the third parameter should be an <intvec>");
    62526151        }
    62536152      }
    62546153      else
    6255       { "ERROR:   wrong list of parameters";
    6256         return();
     6154      { ERROR("ERROR:   wrong list of parameters");
    62576155      }
    62586156    }
    … …  
    62606158  else
    62616159  { if (size(#)>2)
    6262     { "ERROR:   there are too many parameters";
    6263       return();
     6160    { ERROR("ERROR:   there are too many parameters");
    62646161    }
    62656162    int v=0;
    … …  
    62706167                                       // we should get
    62716168  if (ncols(P)<>n)
    6272   { "ERROR:   The first parameter ought to be the matrix of the primary";
    6273     "         invariants."
    6274     return();
     6169  { ERROR("ERROR:   The first parameter ought to be the matrix of the primary"+newline+
     6170    "         invariants.");
    62756171  }
    62766172  if (ncols(REY)<>n)
    6277   { "ERROR:   The second parameter ought to be the Reynolds operator."
    6278     return();
     6173  { ERROR("ERROR:   The second parameter ought to be the Reynolds operator.");
    62796174  }
    62806175  if (v && voice==2)
    … …  
    64266321    { int gen_num=size(#)-1;
    64276322      if (gen_num==0)
    6428       { "ERROR:   There are no generators of the finite matrix group given.";
    6429         return();
     6323      { ERROR("ERROR:   There are no generators of the finite matrix group given.");
    64306324      }
    64316325      int v=#[size(#)];
    64326326      for (i=1;i<=gen_num;i++)
    64336327      { if (typeof(#[i])<>"matrix")
    6434         { "ERROR:   These parameters should be generators of the finite matrix group.";
    6435           return();
     6328        { ERROR("ERROR:   These parameters should be generators of the finite matrix group.");
    64366329        }
    64376330        if ((n<>nrows(#[i])) or (n<>ncols(#[i])))
    6438         { "ERROR:   matrices need to be square and of the same dimensions";
    6439           return();
     6331        { ERROR("ERROR:   matrices need to be square and of the same dimensions");
    64406332        }
    64416333      }
    … …  
    64466338      for (i=1;i<=gen_num;i++)
    64476339      { if (typeof(#[i])<>"matrix")
    6448         { "ERROR:   These parameters should be generators of the finite matrix group.";
    6449           return();
     6340        { ERROR("ERROR:   These parameters should be generators of the finite matrix group.");
    64506341        }
    64516342        if ((n<>nrows(#[i])) or (n<>ncols(#[i])))
    6452         { "ERROR:   matrices need to be square and of the same dimensions";
    6453           return();
     6343        { ERROR("ERROR:   matrices need to be square and of the same dimensions");
    64546344        }
    64556345      }
    … …  
    64576347  }
    64586348  else
    6459   { "ERROR:   There are no generators of the finite matrix group given.";
    6460     return();
     6349  { ERROR("ERROR:   There are no generators of the finite matrix group given.");
    64616350  }
    64626351  if (ncols(P)<>n)
    6463   { "ERROR:   The first parameter ought to be the matrix of the primary";
    6464     "         invariants."
    6465     return();
     6352  { ERROR("ERROR:   The first parameter ought to be the matrix of the primary"+newline+
     6353    "         invariants.");
    64666354  }
    64676355  if (v && voice==2)
    … …  
    66176505"
    66186506{ if (size(#)==0)
    6619   { "ERROR:   There are no generators given.";
    6620     return();
     6507  { ERROR("ERROR:   There are no generators given.");
    66216508  }
    66226509  int ch=char(basering);               // the algorithms depend very much on the
    … …  
    66316518  if (typeof(#[size(#)])=="intvec")
    66326519  { if (size(#[size(#)])<>3)
    6633     { "ERROR:   The <intvec> should have three entries.";
    6634       return();
     6520    { ERROR("ERROR:   The <intvec> should have three entries.");
    66356521    }
    66366522    gen_num=size(#)-1;
    66376523    mol_flag=#[size(#)][1];
    66386524    if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0)))
    6639     { "ERROR:   the second component of <intvec> should be >=0";
    6640       return();
     6525    { ERROR("ERROR:   the second component of <intvec> should be >=0");
    66416526    }
    66426527    int interval=#[size(#)][2];
    … …  
    67506635  if (mol_flag==-1)
    67516636  { if (ch==0)
    6752     { "ERROR:   Characteristic 0 can only be simulated in characteristic p>>0.
    6753 ";
    6754       return();
     6637    { ERROR("ERROR:   Characteristic 0 can only be simulated in characteristic p>>0.");
    67556638    }
    67566639    list L=group_reynolds(#[1..gen_num],v);
    … …  
    67696652  else                                 // the user specified that the
    67706653  { if (ch==0)                         // characteristic divides the group order
    6771     { "ERROR:   The characteristic cannot divide the group order when it is 0.
    6772 ";
    6773       return();
     6654    { ERROR("ERROR:   The characteristic cannot divide the group order when it is 0.");
    67746655    }
    67756656    if (v)
    … …  
    68266707"
    68276708{ if (size(#)<2)
    6828   { "ERROR:   There are too few parameters.";
    6829     return();
     6709  { ERROR("ERROR:   There are too few parameters.");
    68306710  }
    68316711  int ch=char(basering);               // the algorithms depend very much on the
    … …  
    68406720  if (typeof(#[size(#)])=="intvec" && typeof(#[size(#)-1])=="int")
    68416721  { if (size(#[size(#)])<>3)
    6842     { "ERROR:   <intvec> should have three entries.";
    6843       return();
     6722    { ERROR("ERROR:   <intvec> should have three entries.");
    68446723    }
    68456724    gen_num=size(#)-2;
    68466725    mol_flag=#[size(#)][1];
    68476726    if (#[size(#)][2]<0 && (ch==0 or (ch<>0 && mol_flag<>0)))
    6848     { "ERROR:   the second component of <intvec> should be >=0";
    6849       return();
     6727    { ERROR("ERROR:   the second component of <intvec> should be >=0");
    68506728    }
    68516729    int interval=#[size(#)][2];
    … …  
    68536731    int max=#[size(#)-1];
    68546732    if (gen_num==0)
    6855     { "ERROR:   There are no generators of a finite matrix group given.";
    6856       return();
     6733    { ERROR("ERROR:   There are no generators of a finite matrix group given.");
    68576734    }
    68586735  }
    … …  
    68666743    }
    68676744   else
    6868     { "ERROR:   If the two last parameters are not <int> and <intvec>, the last";
    6869       "         parameter should be an <int>.";
    6870       return();
     6745    { ERROR("ERROR:   If the two last parameters are not <int> and <intvec>, the last"+newline+
     6746      "         parameter should be an <int>.");
    68716747    }
    68726748  }
    … …  
    68746750  { if (typeof(#[i])=="matrix")
    68756751    { if (nrows(#[i])<>n or ncols(#[i])<>n)
    6876       { "ERROR:   The number of variables of the base ring needs to be the same";
    6877         "         as the dimension of the square matrices";
    6878         return();
     6752      { ERROR("ERROR:   The number of variables of the base ring needs to be the same"+newline+
     6753        "         as the dimension of the square matrices");
    68796754      }
    68806755    }
    68816756    else
    6882     { "ERROR:   The first parameters should be a list of matrices";
    6883       return();
     6757    { ERROR("ERROR:   The first parameters should be a list of matrices");
    68846758    }
    68856759  }
    … …  
    69856859  if (mol_flag==-1)
    69866860  { if (ch==0)
    6987     { "ERROR:   Characteristic 0 can only be simulated in characteristic p>>0.
    6988 ";
    6989       return();
     6861    { ERROR("ERROR:   Characteristic 0 can only be simulated in characteristic p>>0.");
    69906862    }
    69916863    list L=group_reynolds(#[1..gen_num],v);
    … …  
    70046876  else                                 // the user specified that the
    70056877  { if (ch==0)                         // characteristic divides the group order
    7006     { "ERROR:   The characteristic cannot divide the group order when it is 0.
    7007 ";
    7008       return();
     6878    { ERROR("ERROR:   The characteristic cannot divide the group order when it is 0.");
    70096879    }
    70106880    if (v)
    … …  
    70396909"
    70406910{ if (newring=="")
    7041   { "ERROR:   the second parameter may not be an empty <string>";
    7042     return();
     6911  { ERROR("ERROR:   the second parameter may not be an empty <string>");
    70436912  }
    70446913  if (nrows(F)==1)
    … …  
    70686937  }
    70696938  else
    7070   { "ERROR:   the <matrix> may only have one row";
    7071     return();
     6939  { ERROR("ERROR:   the <matrix> may only have one row");
    70726940  }
    70736941}
    … …  
    71437011    }
    71447012    else
    7145     { "ERROR:    the third parameter may either be empty or a <string>.";
    7146       return();
     7013    { ERROR("ERROR:    the third parameter may either be empty or a <string>.");
    71477014    }
    71487015  }
    … …  
    72277094  }
    72287095  else
    7229   { "ERROR:   the <matrix> may only have one row";
    7230     return();
     7096  { ERROR("ERROR:   the <matrix> may only have one row");
    72317097  }
    72327098}
    … …  
    72677133"
    72687134{ if (newring=="")
    7269   { "ERROR:   the third parameter may not be empty a <string>";
    7270     return();
     7135  { ERROR("ERROR:   the third parameter may not be empty a <string>");
    72717136  }
    72727137  degBound=0;
    … …  
    73167181  }
    73177182  else
    7318   { "ERROR:   the <matrix> may only have one row";
    7319     return();
     7183  { ERROR("ERROR:   the <matrix> may only have one row");
    73207184  }
    73217185}
    … …  
    73627226  }
    73637227  else
    7364   { "ERROR:   the <matrix> may only have one row";
    7365     return();
     7228  { ERROR("ERROR:   the <matrix> may only have one row");
    73667229  }
    73677230}
    … …  
    73757238///////////////////////////////////////////////////////////////////////////////
    73767239
    7377 
    7378 proc secondary_charp_degbound (matrix P, int dmax, list #)
    7379 "USAGE:   secondary_charp_degbound(P,G1,...,Gk[,v]);
    7380 @*          P: a 1xn <matrix> with homogeneous primary invariants, where
    7381                n is the number of variables of the basering;
    7382 @*          G1,...,Gk matrices generating a matrix group for which P are homogeneous
    7383             primary invariants.
    7384 @*          v: an optional <int>;
    7385 ASSUME:   The characteristic of basering is not zero
    7386 RETURN:   secondary invariants of the invariant ring (type <matrix>) and
    7387           irreducible secondary invariants (type <matrix>) of degree at most dmax.
    7388 DISPLAY:  information if v does not equal 0
    7389 THEORY:   A basis of invariants of degrees d<dmax+1 is computed using invariant_basis.
    7390           Among the basis elements and their power products we pick a maximal subset
    7391           that is linearly independent modulo the primary invariants (see paper \"Some
    7392           Algorithms in Invariant Theory of Finite Groups\" by Kemper and
    7393           Steel (1997)) using Groebner basis techniques.
    7394 @*        If dmax is at least the Noether number of the invariant ring, then primary
    7395           invariants and the irreducible secondaries found by secondary_charp_degbound
    7396           are algebra generators of the invariant ring. This may work
    7397           faster than an application of secondary_not_cohen_macaulay.
    7398 SEE ALSO: secondary_charp, secondary_char0, secondary_not_cohen_macaulay
    7399 EXAMPLE:  example secondary_charp_molien; shows an example
    7400 "
    7401 {
    7402  //----------  Internal parameters, whose choice might improve the performance -----
    7403   int pieces = 15;   // For generating reducible secondaries, blocks of #pieces# secondaries
    7404                      // are formed.
    7405                      // If this parameter is small, the memory consumption will decrease.
    7406 
    7407   def br=basering;
    7408 
    7409  //---------------- checking input and setting verbose mode -------------------
    7410   int n=nvars(br);                     // n is the number of variables, as well
    7411                                        // as the size of the matrices, as well
    7412                                        // as the number of primary invariants,
    7413                                        // we should get
    7414   if (char(br)==0)
    7415   { "ERROR:   secondary_charp_degbound should only be used with rings of characteristic p>0.";
    7416     return();
    7417   }
    7418   int i;
    7419   if (size(#)==0)
    7420   { "ERROR: There are no generators given.";
    7421     return();
    7422   }
    7423   if (typeof(#[size(#)])=="int")
    7424   { int v=#[size(#)];
    7425   }
    7426   else
    7427   { int v=0;
    7428   }
    7429   list GEN;
    7430   for (i=1;i<size(#);i++)
    7431   { if ((typeof(#[i])=="matrix") and
    7432         (ncols(#[i])==n) and (nrows(#[i])==n))
    7433     { GEN[i] = #[i];
    7434     }
    7435     else
    7436     { "ERROR: The generators should be square matrices of appropriate size.";
    7437        return();
    7438     }
    7439   }
    7440   if (typeof(#[size(#)])<>"int")
    7441   { if ((typeof(#[size(#)])=="matrix") and
    7442         (ncols(#[size(#)])==n) and (nrows(#[size(#)])==n))
    7443     { GEN[i] = #[size(#)];
    7444     }
    7445     else
    7446     { "ERROR: The last optional parameter should be an integer.";
    7447        return();
    7448     }
    7449   }
    7450 
    7451   if (ncols(P)<>n)
    7452   { "ERROR:   The first parameter ought to be the matrix of the primary";
    7453     "         invariants."
    7454     return();
    7455   }
    7456   if (v && voice==2)
    7457   { "";
    7458   }
    7459   if (v)
    7460   { "  In degree 0 we have: 1";
    7461     "";
    7462   }
    7463  //------------------------ initializing variables ----------------------------
    7464   setring br;
    7465   ideal ProdCand;          // contains products of secondary invariants,
    7466                            // i.e., candidates for reducible sec. inv.
    7467   ideal Mul1,Mul2;
    7468 
    7469   int j, counter, irrcounter, d;
    7470   int m = dmax+1;
    7471   int dgb = degBound;
    7472   degBound = 0;
    7473 //  degBound = dmax;
    7474   option(redSB);
    7475   ideal sP = groebner(ideal(P));
    7476                            // This is the only explicit Groebner basis computation!
    7477   ideal Reductor,SaveRed;  // sP union Reductor is a Groebner basis up to degree i-1
    7478   int SizeSave;
    7479 
    7480   list SSort;         // sec. inv. first sorted by degree and then sorted by the
    7481                       // minimal degree of a non-constant invariant factor.
    7482   list ISSort;        // irr. sec. inv. sorted by degree
    7483 
    7484   poly helpP;
    7485   ideal helpI;
    7486   ideal Indicator;        // will tell us which candidates for sec. inv. we can choose
    7487   ideal ReducedCandidates;
    7488   int helpint;
    7489   int k,k2,k3,minD;
    7490   int ii;
    7491   int saveAttr;
    7492   ideal mon,B,IS;              // IS will contain all irr. sec. inv.
    7493 
    7494   for (i=1;i<m;i++)
    7495   { SSort[i] = list();
    7496     for (k=1;k<=i;k++)
    7497     { SSort[i][k]=ideal();
    7498       attrib(SSort[i][k],"size",0);
    7499     }
    7500   }
    7501  //------------------- generating secondary invariants ------------------------
    7502   for (i=2;i<=m;i++)                   // going through the invariants of degree i-1 -
    7503     {                                  // starting with reducibles
    7504       if (v)
    7505       { "  Looking for Power Products...";
    7506       }
    7507       counter = 0;                      // counts the number of secondaries
    7508       Reductor = ideal(0);
    7509       helpint = 0;
    7510       SaveRed = Reductor;
    7511       SizeSave = 0;
    7512       attrib(Reductor,"isSB",1);
    7513       attrib(SaveRed,"isSB",1);
    7514 
    7515 // We start searching for reducible secondary invariants in degree i-1, i.e., those
    7516 // that are power products of irreducible secondary invariants.
    7517 // It suffices to restrict the search at products of one _irreducible_ sec. inv. (Mul1)
    7518 // with some sec. inv. (Mul2).
    7519 // Moreover, we avoid to consider power products twice since we take a product
    7520 // into account only if the minimal degree of a non-constant invariant factor in "Mul2" is not
    7521 // smaller than the degree of "Mul1".
    7522       for (k=1;k<i-1;k++)
    7523       { if (typeof(ISSort[k])<>"none")
    7524         { Mul1 = ISSort[k];
    7525         }
    7526         else
    7527         { Mul1 = ideal(0);
    7528         }
    7529         if (Mul1[1]<>0)
    7530         { for (minD=k;minD<i-k;minD++)
    7531           { for (k2=1;k2 <= ((attrib(SSort[i-k-1][minD],"size")-1) div pieces)+1; k2++)
    7532             { Mul2=ideal(0);
    7533               if (attrib(SSort[i-k-1][minD],"size")>=k2*pieces)
    7534               { for (k3=1;k3<=pieces;k3++)
    7535                 { Mul2[k3] = SSort[i-k-1][minD][((k2-1)*pieces)+k3];
    7536                 }
    7537               }
    7538               else
    7539               { for (k3=1;k3<=(attrib(SSort[i-k-1][minD],"size") mod pieces);k3++)
    7540                 { Mul2[k3] = SSort[i-k-1][minD][((k2-1)*pieces)+k3];
    7541                 }
    7542               }
    7543 //              if (i<10)
    7544                 ProdCand = simplify(Mul1*Mul2,4);
    7545 //              else { ProdCand = Mul1*Mul2;}
    7546               ReducedCandidates = reduce(ProdCand,sP);
    7547                   // sP union SaveRed union Reductor is a homogeneous Groebner basis
    7548                   // up to degree i-1.
    7549                   // We first reduce by sP (which is fixed, so we can do it once for all),
    7550                   // then by SaveRed resp. by Reductor (which is modified during
    7551                   // the computations).
    7552               Indicator = reduce(ReducedCandidates,SaveRed);
    7553                              // If Indicator[ii]==0 then ReducedCandidates it the reduction
    7554                              // of an invariant that is in the algebra generated by primary
    7555                              // invariants and previously computed secondary invariants.
    7556                              // Otherwise ProdCand[ii] can be taken as secondary invariant.
    7557               if (size(Indicator)<>0)
    7558               { for (ii=1;ii<=ncols(ProdCand);ii++)     // going through all the power products
    7559                 { helpP = Indicator[ii];
    7560                   if (helpP <> 0)
    7561                   { counter++;
    7562                     saveAttr = attrib(SSort[i-1][k],"size")+1;
    7563                     SSort[i-1][k][saveAttr] = ProdCand[ii];
    7564                         // By construction, this is a _reducible_ s.i.
    7565                     attrib(SSort[i-1][k],"size",saveAttr);
    7566                     if (v)
    7567                     { "We found",counter, " secondaries in degree",i-1;
    7568                     }
    7569                     Reductor = ideal(helpP);
    7570                         // Lemma: If G is a homogeneous Groebner basis up to degree i-1 and p is a
    7571                         // homogeneous polynomial of degree i-1 then G union NF(p,G) is
    7572                         // a homogeneous Groebner basis up to degree i-1.
    7573                     attrib(Reductor, "isSB",1);
    7574                         // if Reductor becomes too large, we reduce the whole of Indicator by
    7575                         // it, save it in SaveRed, and work with a smaller Reductor. This turns
    7576                         // out to save a little time.
    7577                     Indicator=reduce(Indicator,Reductor);
    7578                     SizeSave++;
    7579                     SaveRed[SizeSave] = helpP;
    7580                     attrib(SaveRed, "isSB",1);
    7581                   }
    7582                 }
    7583               }
    7584 //              attrib(SaveRed, "isSB",1);
    7585             }
    7586           }
    7587         }
    7588       }
    7589 
    7590       // If there are more sec. inv., then these are irreducible!
    7591       if (v)
    7592       { "Searching irreducible secondary invariants in degree ", i-1;
    7593       }
    7594       B=invariant_basis_modS(i-1,SaveRed,GEN); // B contains enough invariant polynomials of degree d
    7595       irrcounter=0;
    7596       j=0;                         // j goes through all of B -
    7597       // Compare the comments on the computation of reducible sec. inv.!
    7598       ReducedCandidates = reduce(B,sP);
    7599       Indicator = reduce(ReducedCandidates,SaveRed);
    7600       while (size(Indicator)>0)
    7601       { j++;
    7602         helpP = Indicator[j];
    7603         if (helpP <>0)                  // B[j] should be added
    7604         { counter++; irrcounter++;
    7605           IS=IS,B[j];
    7606           saveAttr = attrib(SSort[i-1][i-1],"size")+1;
    7607           SSort[i-1][i-1][saveAttr] = B[j];
    7608           attrib(SSort[i-1][i-1],"size",saveAttr);
    7609           if (typeof(ISSort[i-1]) <> "none")
    7610           { ISSort[i-1][irrcounter] = B[j];
    7611           }
    7612           else
    7613           { ISSort[i-1] = ideal(B[j]);
    7614           }
    7615           if (v)
    7616           { "  We found the irreducible sec. inv. "+string(B[j]);
    7617           }
    7618           Reductor = ideal(helpP);
    7619           attrib(Reductor, "isSB",1);
    7620           Indicator=reduce(Indicator,Reductor);
    7621           SizeSave++;
    7622           SaveRed[SizeSave] = helpP;
    7623           attrib(SaveRed, "isSB",1);
    7624         }
    7625         B[j]=0;
    7626         ReducedCandidates[j]=0;
    7627         Indicator[j]=0;
    7628       }
    7629       if (v)
    7630       { "";
    7631       }
    7632     } // for i
    7633   if (v)
    7634   { "  We're done!";
    7635     "";
    7636   }
    7637   degBound = dgb;
    7638 
    7639   // Prepare return:
    7640   int TotalNumber;
    7641   for (k=1;k<m;k++)
    7642   { for (k2=1;k2<=k;k2++)
    7643     { if (typeof(attrib(SSort[k][k2],"size"))=="int")
    7644       { TotalNumber = TotalNumber + attrib(SSort[k][k2],"size");
    7645       }
    7646     }
    7647   }
    7648   matrix S[1][TotalNumber+1];
    7649   S[1,1]=1;
    7650   j=1;
    7651   for (k=1;k<m;k++)
    7652   { for (k2=1;k2<=k;k2++)
    7653     { if (typeof(attrib(SSort[k][k2],"size"))=="int")
    7654      {for (i=1;i<=attrib(SSort[k][k2],"size");i++)
    7655       { j++;
    7656         S[1,j] = SSort[k][k2][i];
    7657       }
    7658       SSort[k][k2]=ideal();
    7659      }
    7660     }
    7661   }
    7662   return(matrix(S),matrix(compress(IS)));
    7663 }
    7664 example
    7665 { "EXAMPLE: Sturmfels: Algorithms in Invariant Theory 2.3.7 (changed into char 3)"; echo=2;
    7666          ring R=3,(x,y,z),dp;
    7667          matrix A[3][3]=0,1,0,-1,0,0,0,0,-1;
    7668          matrix P=primary_charp_without(A);
    7669          matrix S,IS=secondary_charp_degbound(P,4,A,1);
    7670          print(S);
    7671          print(IS);
    7672 }
    7673 
Note: See TracChangeset for help on using the changeset viewer.