Changeset f86ef8 in git


Ignore:
Timestamp:
Jul 4, 2007, 3:14:43 PM (17 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
2a329dc895a712995f2f237070af97a4eb2ca88a
Parents:
f7e74d328568002cb2b76a79f848a1d7f131710f
Message:
*hannes: p-adic pStd integrated


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

Legend:

Unmodified
Added
Removed
  • Singular/LIB/modstd.lib

    rf7e74d3 rf86ef8  
    11//GP, last modified 23.10.06
    22///////////////////////////////////////////////////////////////////////////////
    3 version="$Id: modstd.lib,v 1.12 2007-06-22 15:30:11 Singular Exp $";
     3version="$Id: modstd.lib,v 1.13 2007-07-04 13:14:43 Singular Exp $";
    44category="Commutative Algebra";
    55info="
     
    88@*       G. Pfister      pfister@mathematik.uni-kl.de
    99@*       H. Schoenemann  hannes@mathematik.uni-kl.de
     10@*       Cindy Magin,    c.magin@web.de
    1011
    1112NOTE:
    1213 A library for computing the Grobner basis of an ideal in the polynomial
    13  ring over the rational numbers using modular methods.The procedures are
     14 ring over the rational numbers using modular methods. The procedures are
    1415 inspired by the following paper:
    1516 Elizabeth A. Arnold:
     
    2324modS(I,L);     liftings to Q of standard bases of I mod p for p in L
    2425primeList(n);  intvec of n primes  <= 2134567879 in decreasing order
     26pStd(p,i);     compute a standard basis of i using p-adic methods
    2527";
    2628
     
    541543}
    542544///////////////////////////////////////////////////////////////////////////////
     545proc pStd(int p,ideal i)
     546"USAGE:  pStd(p,i);p integer, i ideal;
     547RETURN:  an ideal G which is the groebner base for i
     548EXAMPLE: example pStd; shows an example
     549"
     550{
     551  def r=basering;
     552  list rl=ringlist(r);
     553  rl[1]=p;
     554  def r1=ring(rl);
     555  setring r1;
     556  option(redSB);
     557  ideal j=fetch(r,i);
     558  ideal GP=groebner(j);
     559  setring r;
     560  ideal G=fetch(r1,GP);
     561  attrib(G,"isSB",1);
     562  matrix Z=transmat(p,i,G);
     563  matrix G1=gstrich1(p,Z,i,G);
     564  ideal g1=G1;
     565  ideal g22=reduce(g1,G);
     566  matrix G22=transpose(matrix(g22));
     567  matrix M=redmat(G,G1,G22);
     568  matrix Z2=-M*Z;
     569  kill r1;
     570  number c=p;
     571  matrix G0=transpose(matrix(G));
     572  G0= MmodN(G0+ (c)* G22,c^2);
     573  matrix GF=fareyMatrix(G0,c^2);
     574  Z=MmodN(Z+(c)*Z2,c^2);
     575  matrix C=transpose(G);
     576  int n=3;
     577  while(GF<>C)
     578  {
     579    C=GF;
     580    G1= gstrich2(c,Z,i,G0,n);
     581    g1=G1;
     582    g22=reduce(g1,G);
     583    G22=transpose(matrix(g22));
     584    M=redmat(G,G1,G22);
     585    Z2=-M*Z;
     586    Z=MmodN(Z+(c^(n-1))*Z2,c^n);
     587    G0= MmodN(G0+ (c^(n-1))* G22,c^n);
     588    GF=fareyMatrix(G0,c^n);
     589    n++;
     590  }
     591  return(ideal(GF));
     592}
     593example
     594{ "EXAMPLE:"; echo = 2;
     595   ring r=0,(x,y,z),dp;
     596   ideal I=3x3+x2+1,11y5+y3+2,5z4+z2+4;
     597   ideal J=pStd(32003,I);
     598   J;
     599}
     600///////////////////////////////////////////////////////////////////////////
     601proc transmat(int p,ideal i,ideal G)
     602"USAGE:  transmat(p,I,G); p integer, I,G ideal;
     603RETURN:  the transformationmatrix Z for the ideal i mod p and the groebner base for i mod p
     604EXAMPLE: example transmit; shows an example
     605"
     606{
     607  def r=basering;
     608  int n=nvars(r);
     609  list rl=ringlist(r);
     610  rl[1]=p;
     611  def r1=ring(rl);
     612  setring r1;
     613  ideal i=fetch(r,i);
     614  ideal G=fetch(r,G);
     615  attrib(G,"isSB",1);
     616  ring rhelp=p,x(1..n),dp;
     617  list lhelp=ringlist(rhelp);
     618  list l=lhelp[3];
     619  setring r;
     620  rl[3]=l;
     621  def r2=ring(rl);
     622  setring r2;
     623  ideal i=fetch(r,i);
     624  option(redSB);
     625  ideal j=std(i);
     626  matrix T=lift(i,j);
     627  setring r1;
     628  matrix T=fetch(r2,T);
     629  ideal j=fetch(r2,j);
     630  matrix M=lift(j,G);
     631  matrix Z=transpose(T*M);
     632  setring r;
     633  matrix Z=fetch(r1,Z);
     634  return(Z);
     635}
     636example
     637{ "EXAMPLE:"; echo = 2;
     638   ring r=0,(x,y,z),dp;
     639   ideal i=3x3+x2+1,11y5+y3+2,5z4+z2+4;
     640   ideal g=x3-60x2-60, z4-36z2+37, y5+33y3+66;
     641   int p=181;
     642   matrix Z=transmat(p,i,g);
     643   Z;
     644}
     645
     646///////////////////////////////////////////////////////////////////////////
     647proc gstrich1(int p, matrix Z, ideal i, ideal gp)
     648"USAGE:  gstrich1 (p,Z,i,gp); p integer, Z matrix, i,gp ideals;
     649RETURN:  a matrix G such that (Z*F-GP)/p, where F and GP are the matrices of the ideals i and gp
     650"
     651{
     652  matrix F=transpose(matrix(i));
     653  matrix GP=transpose(matrix(gp));
     654  matrix G=(Z*F-GP)/p;
     655  return(G);
     656}
     657///////////////////////////////////////////////////////////////////////////
     658proc gstrich2(number p, matrix Z, ideal i, ideal gp, int n)
     659"USAGE:  gstrich2 (p,Z,i,gp,n); p,n integer, Z matrix, i,gp ideals;
     660RETURN:  a matrix G such that (Z*F-GP)/(p^(n-1)), where F and GP are the matrices of the ideals i and gp
     661"
     662{
     663  matrix F=transpose(matrix(i));
     664  matrix GP=transpose(matrix(gp));
     665  matrix G=(Z*F-GP)/(p^(n-1));
     666  return(G);
     667}
     668///////////////////////////////////////////////////////////////////////////
     669proc redmat(ideal i, matrix h, matrix g)
     670"USAGE:  redmat(i,h,g); i ideal , h,g matrices;
     671RETURN:  a matrix M such that i=M*h+g
     672"
     673{
     674  matrix c=h-g;
     675  ideal f=transpose(c);
     676  matrix N=lift(i,f);
     677  matrix M=transpose(N);
     678  return(M);
     679}
     680///////////////////////////////////////////////////////////////////////////
     681proc fareyMatrix(matrix m,number N)
     682"USAGE:  fareyMatrix(m,y); m matrix, y integer;
     683RETURN:  a matrix k of the matrix m with Farey rational numbers a/b as coefficients
     684EXAMPLE: example fareyMatrix; shows an example
     685"
     686{
     687  ideal I=m;
     688  poly result,p;
     689  int i,j;
     690  number n;
     691  for(i=1;i<=size(I);i++)
     692  {
     693    p=I[i];
     694    result=lead(p);
     695    while(1)
     696    {
     697      if (p==0) {break;}
     698      p=p-lead(p);
     699      n=Farey(N,leadcoef(p));
     700      result=result+n*leadmonom(p);
     701    }
     702    I[i]=result;
     703  }
     704  matrix k=transpose(I);
     705  return(k);
     706}
     707example
     708{"EXAMPLE:"; echo = 2;
     709   ring r=0,(x,y,z),dp;
     710   matrix m[3][1]=x3+682794673x2+682794673,z4+204838402z2+819353608,    y5+186216729y3+372433458;
     711   int p=32003;
     712   matrix b=fareyMatrix(m,p^2);
     713   b;
     714}
     715///////////////////////////////////////////////////////////////////////////
     716proc MmodN(matrix Z,number N)
     717"USAGE:  MmodN(Z,N);Z matrix, N number;
     718RETURN:  the matrix Z mod N
     719EXAMPLE: example MmodN;
     720"
     721{
     722  int i,j,k;
     723  poly m,p;
     724  number c;
     725  for(i=1;i<=nrows(Z);i++)
     726  {
     727    for(j=1;j<=ncols(Z);j++)
     728    {
     729      for(k=1;k<=size(Z[i,j]);k++)
     730      {
     731        m=leadmonom(Z[i,j][k]);
     732        c=leadcoef(Z[i,j][k]) mod N;
     733        p=p+c*m;
     734      }
     735      Z[i,j]=p;
     736      p=0;
     737    }
     738  }
     739  return(Z);
     740}
     741example
     742{ "EXAMPLE:"; echo = 2;
     743   ring r = 0,(x,y,z),dp;
     744   matrix m[3][1]= x3+10668x2+10668, z4-12801z2+12802, y5-8728y3+14547;
     745   number p=32003;
     746   matrix b=MmodN(m,p^2);
     747   b;
     748}
     749///////////////////////////////////////////////////////////////////////////////
    543750/*
    544751ring r=0,(x,y,z),lp;
     
    610817
    611818*/
     819
     820/*
     821ring r=0,(x,y,z),lp;
     822poly s1 = 5x3y2z+3y3x2z+7xy2z2;
     823poly s2 = 3xy2z2+x5+11y2z2;
     824poly s3 = 4xyz+7x3+12y3+1;
     825poly s4 = 3x3-4y3+yz2;
     826ideal i =  s1, s2, s3, s4;
     827
     828ring r=0,(x,y,z),lp;
     829poly s1 = 2xy4z2+x3y2z-x2y3z+2xyz2+7y3+7;
     830poly s2 = 2x2y4z+x2yz2-xy2z2+2x2yz-12x+12y;
     831poly s3 = 2y5z+x2y2z-xy3z-xy3+y4+2y2z;
     832poly s4 = 3xy4z3+x2y2z-xy3z+4y3z2+3xyz3+4z2-x+y;
     833ideal i =  s1, s2, s3, s4;
     834
     835ring r=0,(x,y,z),lp;
     836poly s1 = 8x2y2 + 5xy3 + 3x3z + x2yz;
     837poly s2 = x5 + 2y3z2 + 13y2z3 + 5yz4;
     838poly s3 = 8x3 + 12y3 + xz2 + 3;
     839poly s4 = 7x2y4 + 18xy3z2 +  y3z3;
     840ideal i =  s1, s2, s3, s4;
     841
     842int n = 6;
     843ring r = 0,(x(1..n)),lp;
     844ideal i = cyclic(n);
     845ring s=0,(x(1..n),t),lp;
     846ideal i=imap(r,i);
     847i=homog(i,t);
     848
     849ring r=0,(x(1..4),s),(dp(4),dp);
     850poly s1 =1 + s^2*x(1)*x(3) + s^8*x(2)*x(3) + s^19*x(1)*x(2)*x(4);
     851poly s2 = x(1) + s^8 *x(1)* x(2)* x(3) + s^19* x(2)* x(4);
     852poly s3 = x(2) + s^10*x(3)*x(4) + s^11*x(1)*x(4);
     853poly s4 = x(3) + s^4*x(1)*x(2) + s^19*x(1)*x(3)*x(4) +s^24*x(2)*x(3)*x(4);
     854poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4);
     855ideal i =  s1, s2, s3, s4, s5;
     856
     857ring r=0,(x,y,z),ds;
     858int a =16;
     859int b =15;
     860int c =4;
     861int t =1;
     862poly f =x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3+x^(c-2)*y^c*(y2+t*x)^2;
     863ideal i= jacob(f);
     864
     865ring r=0,(x,y,z),ds;
     866int a =25;
     867int b =25;
     868int c =5;
     869int t =1;
     870poly f =x^a+y^b+z^(3*c)+x^(c+2)*y^(c-1)+x^(c-1)*y^(c-1)*z3+x^(c-2)*y^c*(y2+t*x)^2;
     871ideal i= jacob(f),f;
     872
     873ring r=0,(x,y,z),ds;
     874int a=10;
     875poly f =xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
     876ideal i= jacob(f);
     877
     878ring r=0,(x,y,z),ds;
     879int a =6;
     880int b =8;
     881int c =10;
     882int alpha =5;
     883int beta= 5;
     884int t= 1;
     885poly f =x^a+y^b+z^c+x^alpha*y^(beta-5)+x^(alpha-2)*y^(beta-3)+x^(alpha-3)*y^(beta-4)*z^2+x^(alpha-4)*y^(beta-4)*(y^2+t*x)^2;
     886ideal i= jacob(f);
     887
     888*/
Note: See TracChangeset for help on using the changeset viewer.