Changeset 967543 in git


Ignore:
Timestamp:
Jan 9, 2007, 1:44:33 PM (17 years ago)
Author:
Gerhard Pfister <pfister@…>
Branches:
(u'spielwiese', '17f1d200f27c5bd38f5dfc6e8a0879242279d1d8')
Children:
e19002f9deb88fd1159c864fac5fdee7d9e0a293
Parents:
3d4a5a6feb25d26d6b5353dab0cea5db6d585dab
Message:
die Hauptprozeduren Gert-Martins Wuenschen angepasstCVS: ----------------------------------------------------------------------


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/modstd.lib

    r3d4a5a r967543  
    11//GP, last modified 23.10.06
    22///////////////////////////////////////////////////////////////////////////////
    3 version="$Id: modstd.lib,v 1.1 2007-01-09 10:31:49 Singular Exp $";
     3version="$Id: modstd.lib,v 1.2 2007-01-09 12:44:33 pfister Exp $";
    44category="Commutative Algebra";
    55info="
    … …  
    1818
    1919
     20
    2021PROCEDURES:
    21 modStd(I,1);   compute the reduced Groebner basis of I using modular methods
    22 modS(I,L);     liftings to Q of Groebner bases of I mod p for p in L
     22modStd(I);     compute a standard basis of I using modular methods
     23modS(I,L);     liftings to Q of standard bases of I mod p for p in L
    2324primeList(n);  list of n primes  <= 2134567879 in decreasing order
    2425";
    2526
    26 LIB "crypto.lib";
    27 ///////////////////////////////////////////////////////////////////////////////
    28 proc modStd(ideal I,list #)
    29 "USAGE:  modStd(I,[k]); I ideal (an optional integer k)
    30 RETURN:  if # is empty:
    31          an ideal which is with high probability a reduced Groebner basis of I;
    32          it is not tested whether the result is a Groebner basis and
    33          it is not tested whether the result contains I.
    34          if #[1]=1:
    35          a Groebner basis which contains I if no warning appears;
    36          if I is homogeneous it is a Groebner basis of I
    37 NOTE:    the procedure computes the reduced Groebner basis of I (over the
    38          rational numbers) by using  modular methods. If #[1]=1 and a
    39          warning appears then the result is a Groebner basis with no defined
    40          relation to I; this is a sign that not enough prime numbers have
     27LIB "poly.lib";
     28LIB "krypto.lib";
     29///////////////////////////////////////////////////////////////////////////////
     30proc modStd(ideal I)
     31"USAGE:  modStd(I);
     32RETURN:  a standard basis of I if no warning appears;
     33NOTE:    the procedure computes a standard basis of I (over the
     34         rational numbers) by using  modular methods. If a
     35         warning appears then the result is a standard basis with no defined
     36         relation to I; this is a sign that not enough prime numbers have
    4137         been used. For further experiments see procedure modS.
    4238EXAMPLE: example modStd; shows an example
    4339"
    4440{
     41  int aa=timer;
    4542  def R0=basering;
    4643  list rl=ringlist(R0);
    … …  
    5047  }
    5148  int l,j,k,q;
     49  int en=2134567879;
     50  int an=1000000000;
     51  intvec hi,hl,hc,hpl,hpc;
    5252  list T,TT;
    5353  list L=primeList(5);
    54   L[6]=prime(random(1000000000,2000000000));
     54  L[6]=prime(random(an,en));
    5555  ideal J,cT,lT,K;
    56   ideal I0=I;
    57 
     56  ideal I0=I;
     57  int h=homog(I);
     58  if((!h)&&(ord_test(R0)==0))
     59  {
     60     ERROR("input is not homogeneous and ordering is not local");
     61  }
     62  if(h)
     63  {
     64     execute("ring gn="+string(L[6])+",x(1.."+string(nvars(R0))+"),dp;");
     65     ideal I=fetch(R0,I);
     66     ideal J=std(I);
     67     hi=hilb(J,1);
     68     setring R0;
     69  }
    5870  for (j=1;j<=size(L);j++)
    5971  {
    … …  
    6375    ideal I=fetch(R0,I);
    6476    option(redSB);
    65     ideal I1=groebner(I);
     77    if(h)
     78    {
     79       ideal I1=std(I,hi);
     80    }
     81    else
     82    {
     83      if(ord_test(R0)==-1)
     84      {
     85         ideal I1=std(I);
     86      }
     87      else
     88      {
     89         matrix M;
     90         ideal I1=liftstd(I,M);
     91      }
     92    }
    6693    setring R0;
    6794    T[j]=fetch(oro,I1);
    … …  
    7097  //================= delete unlucky primes ====================
    7198  // unlucky iff the leading ideal is wrong
    72   lT=lead(T[size(T)]);
    73   for (j=1;j<size(T);j++)
    74   {
    75      cT=lead(T[j]);
    76      for(k=1;k<=size(lT);k++)
    77      {
    78         if(lT[k]<cT[k]){lT=cT;break;}
    79      }
    80      if(size(lT)<size(cT)){lT=cT;}
    81   }
    82   j=1;
    83   attrib(lT,"isSB",1);
    84   while(j<=size(T))
    85   {
    86      cT=lead(T[j]);
    87      attrib(cT,"isSB",1);
    88      if((size(reduce(cT,lT))!=0)||(size(reduce(lT,cT))!=0))
    89      {
    90         T=delete(T,j);
    91         L=delete(L,j);
    92         j--;
    93      }
    94      j++;
    95   }
     99  list LL=deleteUnluckyPrimes(T,L);
     100  T=LL[1];
     101  L=LL[2];
     102  lT=LL[3];
    96103  //============ now all leading ideals are the same ============
    97104  for(j=1;j<=ncols(T[1]);j++)
    … …  
    101108      TT[k]=T[k][j];
    102109    }
    103     J[j]=liftPoly(TT,L);
     110    J[j]=liftPoly(TT,L); 
    104111  }
    105112  //=========== chooses more primes up to the moment the result becomes stable
    … …  
    107114  {
    108115     k=0;
    109      q=prime(random(2000000011,2100000000));
     116     q=prime(random(an,en));
    110117     while(k<size(L))
    111118     {
    … …  
    114121        {
    115122           k=0;
    116            q=prime(random(1000000,2100000000));
     123           q=prime(random(an,en));
    117124        }
    118125     }
    … …  
    123130     ideal i=fetch(R0,I);
    124131     option(redSB);
    125      i=groebner(i);
     132     if(h)
     133     {
     134       i=std(i,hi);
     135     }
     136     else
     137     {
     138       if(ord_test(R0)==-1)
     139       {
     140          i=std(i);
     141       }
     142       else
     143       {
     144          matrix M;
     145          i=liftstd(i,M);
     146       }
     147     }
    126148     setring R0;
    127149     T[size(T)+1]=fetch(@r,i);
    … …  
    146168        }
    147169        k=1;
    148         for(j=1;j<=size(K);j++){if(K[j]-J[j]!=0){k=0;break;}}
     170        for(j=1;j<=size(K);j++)
     171        {
     172           if(K[j]-J[j]!=0)
     173           {
     174              k=0;
     175              J=K;
     176              break;
     177           }
     178        }
    149179     }
    150180     if(k){break;}
    151181  }
    152   //============ optional test for standard basis and I=J =======
    153   if(size(#)>0)
    154   {
    155     J=std(J);
    156     I0=reduce(I0,J);
    157     if(size(I0)>0){"WARNING: The input ideal is not contained
     182  //============  test for standard basis and I=J =======
     183"computed";timer-aa;aa=timer;
     184  J=std(J);
     185  I0=reduce(I0,J);
     186  if(size(I0)>0){"WARNING: The input ideal is not contained
    158187                        in the ideal generated by the standardbasis";}
    159   }
    160188  attrib(J,"isSB",1);
     189"verification";timer-aa;
    161190  return(J);
    162191}
    … …  
    167196   ideal J=modStd(I);
    168197   J;
    169 }
    170 ///////////////////////////////////////////////////////////////////////////////
    171 proc modS(ideal I, list L)
     198}   
     199///////////////////////////////////////////////////////////////////////////////
     200proc modS(ideal I, list L, list #)
    172201"USAGE:  modS(I,L); I ideal, L list of primes
     202         if size(#)>0 std is used instead of groebner
    173203RETURN:  an ideal which is with high probability a standard basis
    174204NOTE:    This procedure is designed for fast experiments.
    … …  
    195225    ideal i=fetch(R0,I);
    196226    option(redSB);
    197     i=groebner(i);
     227    if(size(#)>0)
     228    {
     229       i=std(i);
     230    }
     231    else
     232    {
     233       i=groebner(i)
     234    }
    198235    setring R0;
    199236    T[j]=fetch(@r,i);
    … …  
    202239  //================= delete unlucky primes ====================
    203240  // unlucky iff the leading ideal is wrong
    204   lT=lead(T[1]);
    205   for (j=2;j<=size(T);j++)
    206   {
    207      cT=lead(T[j]);
    208      for(k=1;k<=size(lT);k++)
    209      {
    210         if(lT[k]<cT[k]){lT=cT;break;}
    211      }
    212      if(size(lT)<size(cT)){lT=cT;}
    213   }
    214   j=1;
    215   attrib(lT,"isSB",1);
    216   while(j<=size(T))
    217   {
    218      cT=lead(T[j]);
    219      attrib(cT,"isSB",1);
    220      if((size(reduce(cT,lT))!=0)||(size(reduce(lT,cT))!=0))
    221      {
    222         T=delete(T,j);
    223         L=delete(L,j);
    224         j--;
    225      }
    226      j++;
    227   }
     241  list LL=deleteUnluckyPrimes(T,L);
     242  T=LL[1];
     243  L=LL[2];
    228244  //============ now all leading ideals are the same ============
    229245  for(j=1;j<=ncols(T[1]);j++)
    … …  
    233249      TT[k]=T[k][j];
    234250    }
    235     J[j]=liftPoly(TT,L);
    236 
     251    J[j]=liftPoly(TT,L); 
    237252  }
    238253  attrib(J,"isSB",1);
    … …  
    246261   ideal J=modS(I,L);
    247262   J;
    248 }
     263} 
     264///////////////////////////////////////////////////////////////////////////////
     265proc deleteUnluckyPrimes(list T,list L)
     266"USAGE:  deleteUnluckyPrimes(T,L);T list of polys, L list of primes 
     267RETURN:  list L,T with T list of polys, L list of primes
     268EXAMPLE: example deleteUnluckyPrimes; shows an example
     269NOTE:    works only for homogeneous ideals with global orderings or
     270         for ideals with local orderings
     271"
     272{
     273  int j,k;
     274  intvec hl,hc;
     275  ideal cT,lT;
     276
     277  lT=lead(T[size(T)]);
     278  attrib(lT,"isSB",1);
     279  hl=hilb(lT,1);
     280  for (j=1;j<size(T);j++)
     281  {
     282     cT=lead(T[j]);
     283     attrib(cT,"isSB",1);
     284     hc=hilb(cT,1);
     285     if(hl==hc)
     286     {
     287        for(k=1;k<=size(lT);k++)
     288        {
     289           if(lT[k]<cT[k]){lT=cT;break;}
     290           if(lT[k]>cT[k]){break;}
     291        }
     292     }
     293     else
     294     {
     295        if(hc<hl){lT=cT;hl=hilb(lT,1);}
     296     }
     297  }
     298  j=1;
     299  attrib(lT,"isSB",1);
     300  while(j<=size(T))
     301  {
     302     cT=lead(T[j]);
     303     attrib(cT,"isSB",1);
     304     if((size(reduce(cT,lT))!=0)||(size(reduce(lT,cT))!=0))
     305     {
     306        T=delete(T,j);
     307        L=delete(L,j);
     308        j--;
     309     }
     310     j++;
     311  }
     312  return(list(T,L,lT));
     313}
     314example
     315{ "EXAMPLE:"; echo = 2;
     316   list L=2,3,5,7,11;
     317   ring r=0,(y,x),Dp;
     318   ideal I1=y2x,y6;
     319   ideal I2=yx2,y3x,x5,y6;
     320   ideal I3=y2x,x3y,x5,y6;
     321   ideal I4=y2x,x3y,x5;
     322   ideal I5=y2x,yx3,x5,y6;
     323   list T=I1,I2,I3,I4,I5;
     324   list TT=deleteUnluckyPrimes(T,L);
     325   TT;
     326} 
    249327///////////////////////////////////////////////////////////////////////////////
    250328proc liftPoly(list T, list L)
    … …  
    259337   list TT;
    260338   number n;
    261 
     339 
    262340   number N=L[1];
    263341   for(i=2;i<=size(L);i++)
    … …  
    288366       }
    289367     }
    290      n=chineseR(TT,L);
    291      n=Farey(N,n);
     368     n=chineseR(TT,L,N);
     369     n=Farey(N,n); 
    292370     result=result+n*p;
    293371   }
    … …  
    300378   liftPoly(T,L);
    301379}
     380 ///////////////////////////////////////////////////////////////////////////////
     381proc fareyIdeal(ideal I,list L)
     382{
     383   poly result,p;
     384   int i,j;
     385   number n;
     386   number N=L[1];
     387   for(i=2;i<=size(L);i++)
     388   {
     389      N=N*L[i];
     390   }
     391
     392   for(i=1;i<=size(I);i++)
     393   {
     394     p=I[i];
     395     result=lead(p);
     396     while(1)
     397     {
     398        if (p==0) {break;}
     399        p=p-lead(p);
     400        n=Farey(N,leadcoef(p));
     401        result=result+n*leadmonom(p);
     402     }
     403     I[i]=result;
     404   }
     405   return(I);
     406}
    302407///////////////////////////////////////////////////////////////////////////////
    303408proc Farey (number P, number N)
    304409"USAGE:  Farey (P,N); P, N number;
    305 RETURN:  a rational number a/b such that a/b=N mod P
     410RETURN:  a rational number a/b such that a/b=N mod P 
    306411         and |a|,|b|<(P/2)^{1/2}
    307412"
    … …  
    314419   while (N!=0)
    315420   {
    316         if (2*N^2<P){return(N/B);}
    317         D=E mod N;
    318         C=A-(E-E mod N)/N*B;
    319         E=N;
    320         N=D;
    321         A=B;
    322         B=C;
     421        if (2*N^2<P)
     422        {           
     423           return(N/B);
     424        }
     425        D=E mod N;
     426        C=A-(E-E mod N)/N*B;
     427        E=N;
     428        N=D;
     429        A=B;
     430        B=C;
    323431   }
    324432   return(0);
    … …  
    329437   Farey(32003,12345);
    330438}
    331 ///////////////////////////////////////////////////////////////////////////////
    332 proc chineseR(list T,list L)
    333 "USAGE:  chineseRem(T,L);
     439 ///////////////////////////////////////////////////////////////////////////////
     440proc chineseR(list T,list L,number N)
     441"USAGE:  chineseR(T,L);
    334442RETURN: x such that x = T[i] mod L[i]
    335443NOTE:   chinese remainder theorem
    336 EXAMPLE:example chineseRem; shows an example
     444EXAMPLE:example chineseR; shows an example
    337445"
    338446{
    … …  
    341449   {
    342450      x=T[1] mod L[1];
    343       if(x>L[1]/2){x=x-L[1];}
    344451      return(x);
    345452   }
    346453   int i;
    347454   int n=size(L);
    348    number N=1;
    349    for(i=1;i<=n;i++)
    350    {
    351       N=N*L[i];
    352    }
    353455   list M;
    354456   for(i=1;i<=n;i++)
    … …  
    362464   }
    363465   x=x mod N;
    364    if (x>N/2) { x=x-N; }
    365466   return(x);
    366467}
    … …  
    368469{ "EXAMPLE:"; echo = 2;
    369470   ring R = 0,x,dp;
    370    chineseRem(list(24,15,7),list(2,3,5));
    371 }
     471   chineseR(list(24,15,7),list(2,3,5));
     472}
     473 
    372474///////////////////////////////////////////////////////////////////////////////
    373475proc primeList(int n)
    374476"USAGE:  primeList(n);
    375 RETURN: the list of n greatest primes  <= 2134567879
     477RETURN: the list of n greatest primes  <= 2134567879   
    376478EXAMPLE:example primList; shows an example
    377479"
    … …  
    392494   size(L);
    393495   L[size(L)];
    394 }
     496}   
    395497///////////////////////////////////////////////////////////////////////////////
    396498/*
    … …  
    416518ideal i =  s1, s2, s3, s4;
    417519
    418 ring r=0,x(1..4),lp;
    419 ideal i=cyclic(4);
     520int n = 6;
     521ring r = 0,(x(1..n)),lp;
     522ideal i = cyclic(n);
     523ring s=0,(x(1..n),t),lp;
     524ideal i=imap(r,i);
     525i=homog(i,t);
    420526
    421527ring r=0,(x(1..4),s),(dp(4),dp);
    … …  
    426532poly s5 = x(4) + s^31* x(1)* x(2)* x(3)* x(4);
    427533ideal i =  s1, s2, s3, s4, s5;
     534
     535ring r=0,(x,y,z),ds;
     536int a =16;
     537int b =15;
     538int c =4;
     539int t =1;
     540poly 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;
     541ideal i= jacob(f);
     542
     543ring r=0,(x,y,z),ds;
     544int a =25;
     545int b =25;
     546int c =5;
     547int t =1;
     548poly 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;
     549ideal i= jacob(f),f;
     550
     551ring r=0,(x,y,z),ds;
     552int a=10;
     553poly f =xyz*(x+y+z)^2 +(x+y+z)^3 +x^a+y^a+z^a;
     554ideal i= jacob(f);
     555
     556ring r=0,(x,y,z),ds;
     557int a =6;
     558int b =8;
     559int c =10;
     560int alpha =5;
     561int beta= 5;
     562int t= 1;
     563poly 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;
     564ideal i= jacob(f);
     565
    428566*/
Note: See TracChangeset for help on using the changeset viewer.