Changeset 3517649 in git


Ignore:
Timestamp:
Jun 12, 2009, 10:01:11 AM (15 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
008e8143a13845539fcd39e47a9b1267e6d78265
Parents:
b5d4d19ebfd46b641a2428d8d8f86d9d6dafa6da
Message:
*pfister: fixes


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/modstd.lib

    rb5d4d1 r3517649  
    11//GP, last modified 23.10.06
    22///////////////////////////////////////////////////////////////////////////////
    3 version="$Id: modstd.lib,v 1.15 2009-04-15 11:14:36 seelisch Exp $";
     3version="$Id: modstd.lib,v 1.16 2009-06-12 08:01:11 Singular Exp $";
    44category="Commutative Algebra";
    55info="
    … …  
    2828LIB "crypto.lib";
    2929///////////////////////////////////////////////////////////////////////////////
     30proc pTestSB(ideal I, ideal J, list L)
     31"USAGE:  pTestSB(I,J,L); I,J ideals, L intvec of primes
     32RETURN: 1 (resp. 0) if for a randomly chosen prime p not in L
     33        J mod p is (resp. is not) a standard basis of I mod p
     34EXAMPLE:example primList; shows an example
     35"
     36{
     37  int i,j,k,p;
     38  def R=basering;
     39  list r= ringlist(R);
     40
     41  while(!j)
     42  {
     43    j=1;
     44    p=prime(random(1000000000,2134567879));
     45    for(i=1;i<=size(L);i++)
     46    {
     47      if(p==L[i]){j=0;break}
     48    }
     49    if(j)
     50    {
     51      for(i=1;i<=ncols(J);i++)
     52      {
     53        for(k=2;k<=size(J[i]);k++)
     54        {
     55          if((denominator(leadcoef(J[i][k])) mod
     56p)==0){j=0;break}
     57        }
     58        if(!j){break;}
     59      }
     60    }
     61  }
     62  r[1]=p;
     63  def @R=ring(r);
     64  setring @R;
     65  ideal I=imap(R,I);
     66  ideal J=imap(R,J);
     67  attrib(J,"isSB",1);
     68  j=1;
     69  if(size(reduce(I,J))!=0){j=0;}
     70  if(j)
     71  {
     72    ideal K=std(J);
     73    if(size(reduce(K,J))!=0){j=0;}
     74  }
     75  setring R;
     76  return(j);
     77}
     78example
     79{ "EXAMPLE:"; echo = 2;
     80   intvec L=2,3,5;
     81   ring r=0,(x,y,z),dp;
     82   ideal i=x+1,x+y+1;
     83   ideal j=x+1,y;
     84   pTestSB(i,i,L);
     85   pTestSB(i,j,L);
     86}
     87
     88proc primeList(int n, list #)
     89"USAGE:  primeList(n); (resp. primeList(n,L); )
     90RETURN: the intvec of n greatest primes  <= 2147483647 (resp. n greatest
     91        primes < L[size(L)] union with L)
     92EXAMPLE:example primList; shows an example
     93"
     94{
     95  intvec L;
     96  int i,p;
     97  if(size(#)==0)
     98  {
     99      p=2147483647;
     100      L[1]=p;
     101   }
     102   else
     103  {
     104     L=#[1];
     105     p=prime(L[size(L)]-1);
     106     L[size(L)+1]=p;
     107  }
     108  if(p==2){ERROR("no more primes");}
     109  for(i=2;i<=n;i++)
     110  {
     111    p=prime(p-1);
     112    L[size(L)+1]=p;
     113  }
     114  return(L);
     115}
     116example
     117{ "EXAMPLE:"; echo = 2;
     118   intvec L=primeList(10);
     119   size(L);
     120   L[size(L)];
     121   L=primeList(5,L);
     122   size(L);
     123   L[size(L)];
     124}
     125
     126
    30127proc modStd(ideal I)
    31128"USAGE:  modStd(I);
    … …  
    33130NOTE:    the procedure computes a standard basis of I (over the
    34131         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
    37          been used. For further experiments see procedure modS.
     132         warning appears then the result is a standard basis
     133         containing I and with high probability a standard basis of I.
     134         For further experiments see procedure modS.
    38135EXAMPLE: example modStd; shows an example
    39136"
    … …  
    43140  if((npars(R0)>0)||(rl[1]>0))
    44141  {
    45      ERROR("characteristic of basering should be zero");
    46   }
    47   int l,j,k,q;
     142     ERROR("characteristic of basering should be zero, basering should have
     143no parameters");
     144  }
     145
     146  int k,c;
     147  int pd=printlevel-voice+2;
     148  int j=1;
     149  int h=homog(I);
    48150  int en=2134567879;
    49151  int an=1000000000;
    50   intvec hi,hl,hc,hpl,hpc;
    51   list T,TT;
    52   intvec L=primeList(5);
    53   L[6]=prime(random(an,en));
    54   ideal J,cT,lT,K;
    55   ideal I0=I;
    56   int h=homog(I);
    57   if((!h)&&(ord_test(R0)==0))
    58   {
    59      ERROR("input is not homogeneous and ordering is not local");
    60   }
    61   if(h)
    62   {
    63      execute("ring gn="+string(L[6])+",x(1.."+string(nvars(R0))+"),dp;");
    64      ideal I=fetch(R0,I);
    65      ideal J=std(I);
    66      hi=hilb(J,1);
    67      setring R0;
    68   }
    69   for (j=1;j<=size(L);j++)
    70   {
    71     rl[1]=L[j];
    72     def oro=ring(rl);
    73     setring oro;
    74     ideal I=fetch(R0,I);
    75     option(redSB);
    76     if(h)
    77     {
    78        ideal I1=std(I,hi);
    79     }
    80     else
    81     {
    82       if(ord_test(R0)==-1)
    83       {
    84          ideal I1=std(I);
    85       }
    86       else
    87       {
    88          matrix M;
    89          ideal I1=liftstd(I,M);
    90       }
    91     }
    92     setring R0;
    93     T[j]=fetch(oro,I1);
    94     kill oro;
    95   }
     152
     153  intvec opt = option(get);     // Save current options
     154  intvec L=primeList(10);
     155  L[5]=prime(random(an,en));
     156  list T,TT,TH,LL;
     157
     158
     159  ideal J,K;
     160
     161  option(redSB);
     162
     163  while(1)
     164  {
     165     while(j<=size(L))
     166     {
     167        if(pd>2){c++;c;}
     168        rl[1]=L[j];
     169        def @r=ring(rl);
     170        setring @r;
     171        ideal i=fetch(R0,I);
     172        i=groebner(i);
     173        setring R0;
     174        T[size(T)+1]=fetch(@r,i);
     175        kill @r;
     176        j++;
     177     }
     178     if(pd>2){"lifting";}
    96179  //================= delete unlucky primes ====================
    97180  // unlucky if and only if the leading ideal is wrong
    98   list LL=deleteUnluckyPrimes(T,L);
    99   T=LL[1];
    100   L=LL[2];
    101   lT=LL[3];
     181     LL=deleteUnluckyPrimes(T,L);
     182     TH=LL[1];
     183     L=LL[2];
    102184  //============ now all leading ideals are the same ============
    103   for(j=1;j<=ncols(T[1]);j++)
    104   {
    105     for(k=1;k<=size(L);k++)
    106     {
    107       TT[k]=T[k][j];
    108     }
    109     J[j]=liftPoly(TT,L);
    110   }
    111   //=========== chooses more primes up to the moment the result becomes stable
    112   while(1)
    113   {
    114      k=0;
    115      q=prime(random(an,en));
    116      while(k<size(L))
     185     for(j=1;j<=ncols(TH[1]);j++)
    117186     {
    118         k++;
    119         if(L[k]==q)
     187         for(k=1;k<=size(L);k++)
     188         {
     189             TT[k]=TH[k][j];
     190         }
     191         J[j]=liftPoly(TT,L);
     192     }
     193     if(pd>2){"list of primes:";L;"pTest";}
     194     if(pTestSB(I,J,L))
     195     {
     196        attrib(J,"isSB",1);
     197        K=std(J);
     198        if(size(reduce(I,J))==0)
    120199        {
    121            k=0;
    122            q=prime(random(an,en));
    123         }
    124      }
    125      L[size(L)+1]=q;
    126      rl[1]=L[size(L)];
    127      def @r=ring(rl);
    128      setring @r;
    129      ideal i=fetch(R0,I);
    130      option(redSB);
    131      if(h)
    132      {
    133        i=std(i,hi);
    134      }
    135      else
    136      {
    137        if(ord_test(R0)==-1)
    138        {
    139           i=std(i);
    140        }
    141        else
    142        {
    143           matrix M;
    144           i=liftstd(i,M);
    145        }
    146      }
    147      setring R0;
    148      T[size(T)+1]=fetch(@r,i);
    149      kill @r;
    150      cT=lead(T[size(T)]);
    151      attrib(cT,"isSB",1);
    152      if((size(reduce(cT,lT))!=0)||(size(reduce(lT,cT))!=0))
    153      {
    154         T=delete(T,size(T));
    155         L=L[1..size(L)-1];
    156         k=0;
    157      }
    158      else
    159      {
    160         for(j=1;j<=ncols(T[1]);j++)
    161         {
    162            for(k=1;k<=size(L);k++)
     200           if(size(reduce(K,J))==0)
    163201           {
    164               TT[k]=T[k][j];
     202               if(!((h)||(ord_test(R0)==-1)))
     203              {
     204
     205"==================================================================";
     206                 "    The input is not homogeneous and the ordering is not
     207local.   ";
     208                 "WARNING: ideal generated by output may be greater then
     209input ideal";
     210
     211"==================================================================";
     212              }
     213              option(set, opt);
     214              return(J);
    165215           }
    166            K[j]=liftPoly(TT,L);
    167         }
    168         k=1;
    169         for(j=1;j<=size(K);j++)
    170         {
    171            if(K[j]-J[j]!=0)
    172            {
    173               k=0;
    174               J=K;
    175               break;
    176            }
    177         }
    178      }
    179      if(k){break;}
    180   }
    181   //============  test for standard basis and I=J =======
    182   J=std(J);
    183   I0=reduce(I0,J);
    184   if(size(I0)>0)
    185   {
    186      "WARNING: The input ideal is not contained
    187                         in the ideal generated by the standardbasis";
    188      "list of primes used:";
    189      L;
    190   }
    191   attrib(J,"isSB",1);
    192   return(J);
    193 }
    194 example
    195 { "EXAMPLE:"; echo = 2;
    196    ring r=0,(x,y,z),dp;
     216           if(pd>2){"pTest o.k. but result wrong";}
     217         }
     218         if(pd>2){"pTest o.k. but result wrong";}
     219      }
     220      j=size(L)+1;
     221      L=primeList(10,L);
     222   }
     223}
     224example
     225{ "EXAMPLE:"; echo = 2;
     226   ring r=0,(x,y,z,t),dp;
    197227   ideal I=3x3+x2+1,11y5+y3+2,5z4+z2+4;
     228   ideal J=modStd(I);
     229   J;
     230   I=homog(I,t);
     231   J=modStd(I);
     232   J;
     233   ring s=0,(x,y,z),ds;
     234   ideal I=jacob(x5+y6+z7+xyz);
    198235   ideal J=modStd(I);
    199236   J;
    … …  
    517554}
    518555
    519 ///////////////////////////////////////////////////////////////////////////////
    520 proc primeList(int n)
    521 "USAGE:  primeList(n);
    522 RETURN: the intvec of n greatest primes  <= 2134567879
    523 EXAMPLE:example primList; shows an example
    524 "
    525 {
    526   intvec L=0:n;
    527   int i;
    528   int p=2134567879;
    529   for(i=1;i<=n;i++)
    530   {
    531     L[i]=p;
    532     p=prime(p-1);
    533   }
    534   return(L);
    535 }
    536 example
    537 { "EXAMPLE:"; echo = 2;
    538    intvec L=primeList(10);
    539    size(L);
    540    L[size(L)];
    541 }
    542556///////////////////////////////////////////////////////////////////////////////
    543557proc pStd(int p,ideal i)
Note: See TracChangeset for help on using the changeset viewer.