Changeset 8b87364 in git

Timestamp:

Nov 5, 2001, 5:06:29 PM (22 years ago)

Author:

Gerhard Pfister <pfister@…>

Branches:

(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')

Children:

16105c8f405ac4e9d787fbd2bae347a9686645db

Parents:

5b3302cc6f7c9964ee19073aca0b273ad9e606a9

Message:

neue proceduren eingefuegt


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

Location:

Singular/LIB

Files:

• general.lib (modified) (3 diffs)
• presolve.lib (modified) (3 diffs)
• standard.lib (modified) (3 diffs)

Singular/LIB/general.lib

@@ -2 +2 @@
 //anne, added deleteSublist and watchdog 12.12.2000
 ///////////////////////////////////////////////////////////////////////////////
-version="$Id: general.lib,v 1.40 2001-08-27 14:47:50 Singular Exp $";
+version="$Id: general.lib,v 1.41 2001-11-05 16:05:48 pfister Exp $";
 category="General purpose";
 info="
@@ -24 +24 @@
  watchdog(i,cmd);       only wait for result of command cmd for i seconds
  which(command);        search for command and return absolute path, if found
+ primecoeffs(J[,q]);    primefactors <= min(p,32003) of coeffs of J 
+ primefactors(n [,p]);  primefactors <= min(p,32003) of n
            (parameters in square brackets [] are optional)
 ";
 
 LIB "inout.lib";
+LIB "poly.lib";
+LIB "matrix.lib";
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -1211 +1215 @@
 }
 ///////////////////////////////////////////////////////////////////////////////
+proc primefactors (n, list #)
+"USAGE:   primefactors(n [,p]); n = int or number, p = integer
+COMPUTE: primefactors <= min(p,32003) of n (default p = 32003)
+RETURN:  a list, say l, 
+         l[1] : primefactors <= min(p,32003) of n 
+         l[2] : l[2][i] = multiplicity of l[1][i]
+         l[3] : remaining factor ( n=product{ (l[1][i]^l[2][i])*l[3]} )
+                type(l[3])=typeof(n)
+NOTE:    If n is a long integer (of type number) then the procedure
+         finds primefactors <= min(p,32003) but n may be larger as
+         2147483647 (max. integer representation)
+WARNING: the procedure works for small integers only, just by testing all
+         primes (not to be considerd as serious prime factorization!)
+EXAMPLE: example primefactors; shows an example
+"
+{
+   int ii,jj,z,p,num,w3,q;
+   intvec w1,w2,v;
+   list l;
+   if (size(#) == 0) 
+   {
+      p=32003; 
+   }
+   else 
+   {
+     if( typeof(#[1]) != "int")
+     {
+       ERROR("2nd parameter must be of type int"+newline);
+     }
+     p=#[1];
+   }
+   if( n<0) { n=-n;};
+
+// ----------------- case: 1st parameter is a number --------------------   
+   if (typeof(n) =="number")
+   {
+     kill w3;
+     number w3;
+     if( n > 2147483647 )           //2147483647 max. integer representation
+     {
+        v = primes(2,p);
+        number m;
+        for( ii=1; ii<=size(v); ii++)
+        { 
+          jj=0;
+          while(1)
+          { 
+            q  = v[ii];
+            jj = jj+1;             
+            m  = n/q;                  //divide n as often as possible
+            if (denominator(m)!=1) { break;  }
+            n=m;
+          }
+          if( jj>1 ) 
+          {
+             w1 = w1,v[ii];          //primes
+             w2 = w2,jj-1;           //powers
+          }
+          if( n <= 2147483647 ) { break;  }
+        }
+     }
+
+     if( n >  2147483647 )            //n is still too big
+     {
+       if( size(w1) >1 )               //at least 1 primefactor was found
+       {
+         w1 = w1[2..size(w1)];
+         w2 = w2[2..size(w2)];
+       } 
+       else                           //no primefactor was found
+       {
+         w1 = 1; w2 = 1;
+       }     
+       l  = w1,w2,n;
+       return(l);
+     }
+
+     if( n <= 2147483647 )          //n is in inter range
+     {
+       num = int(n);
+       kill n;
+       int n = num;
+     }
+   }
+   
+// --------------------------- trivial cases --------------------
+   if( n==0 ) 
+   {  
+     w1=1; w2=1; w3=0; l=w1,w2,w3;
+     return(l);
+   }
+   
+   if( n==1 ) 
+   {   
+       w3=1;
+       if( size(w1) >1 )               //at least 1 primefactor was found
+       {
+         w1 = w1[2..size(w1)];
+         w2 = w2[2..size(w2)];
+       } 
+       else                           //no primefactor was found
+       {
+         w1 = 1; w2 = 1;
+       }     
+      l=w1,w2,w3;
+      return(l);
+   }
+   if ( prime(n)==n )         //note: prime(n) <= 32003 in Singular
+   {                          //case n is a prime
+      if (p > n)
+      {  
+        w1=w1,n; w2=w2,1; w3=1;
+        w1 = w1[2..size(w1)];
+        w2 = w2[2..size(w2)];
+        l=w1,w2,w3;
+        return(l);
+      }
+      else
+      {
+        w3=n;
+        if( size(w1) >1 )               //at least 1 primefactor was found
+        {
+           w1 = w1[2..size(w1)];
+           w2 = w2[2..size(w2)];
+         } 
+         else                           //no primefactor was found
+         {
+           w1 = 1; w2 = 1;
+         }     
+         l=w1,w2,w3;
+        return(l);
+      }      
+   }
+   else 
+   {
+      if ( p >= n)
+      {
+         v = primes(q,n div 2 + 1);
+      }
+      else
+      {
+         v = primes(q,p);
+      }
+//------------- search for primfactors <= last entry of v ------------    
+      for(ii=1; ii<=size(v); ii++)
+      {
+         z=0;
+         while( (n mod v[ii]) == 0 )
+         {  
+            z=z+1;
+            n = n div v[ii];
+         }
+         if (z!=0)
+         { 
+            w1 = w1,v[ii];        //primes
+            w2 = w2,z;            //multiplicities
+         }
+      }
+   }
+//--------------- case:at least 1 primefactor was found ---------------
+   if( size(w1) >1 )               //at least 1 primefactor was found
+   {
+      w1 = w1[2..size(w1)];
+      w2 = w2[2..size(w2)];
+   } 
+   else                           //no primefactor was found
+   {
+     w1 = 1; w2 = 1;
+   } 
+   w3 = n;
+   l  = w1,w2,w3;
+   return(l);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   primefactors(7*8*121);
+   ring r = 0,x,dp;
+   primefactors(123456789100);
+}  
+
+///////////////////////////////////////////////////////////////////////////////
+proc primecoeffs(J, list #)
+"USAGE:   primecoeffs(J[,q]); J any type which can be converted to a matrix
+         e.g. ideal, matrix, vector, module, int, intvec
+         q = intger
+COMPUTE: primefactors <= min(p,32003) of coeffs of J (default p = 32003)
+RETURN:  a list, say l, of two intvectors:
+         l[1] : the different primefactors of all coefficients of J
+         l[2] : the different remaining factors
+NOTE:    the procedure works for small integers only, just by testing all
+         primes (not to be considerd as serious prime factorization!)
+EXAMPLE: example primecoeffs; shows an example
+"
+{
+   int q,ii,n,mark;;
+   if (size(#) == 0) 
+   {
+      q=32003; 
+   }
+   else 
+   {
+     if( typeof(#[1]) != "int")
+     {
+       ERROR("2nd parameter must be of type int"+newline);
+     }
+     q=#[1];
+   }
+
+   if (defined(basering) == 0)
+   {
+     mark=1;
+     ring r = 0,x,dp;
+   }
+   def I = ideal(matrix(J));
+   poly p = product(maxideal(1));
+   matrix Coef=coef(I[1],p); 
+   ideal id, jd, rest;
+   intvec v,re;
+   list result,l;
+   for(ii=2; ii<=ncols(I); ii++)
+   {
+     Coef=concat(Coef,coef(I[ii],p));
+   }
+   id = Coef[2,1..ncols(Coef)];
+   id = simplify(id,6);
+   for (ii=1; ii<=size(id); ii++)  
+   {   
+     l = primefactors(number(id[ii]),q);     
+     jd = jd,l[1];
+     rest = rest,l[3];
+   } 
+   jd = simplify(jd,6);
+   for (ii=1; ii<=size(jd); ii++)  
+   { 
+     v[ii]=int(jd[ii]);
+   }
+   v = sort(v)[1];
+   rest = simplify(rest,6);
+   id = sort(id)[1];
+   if (mark)
+   {
+     for (ii=1; ii<=size(rest); ii++)
+     {
+       re[ii] = int(rest[ii]);
+     }
+     result = v,re;
+   }
+   else
+   {
+      result = v,rest; 
+   }
+   return(result);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   primecoeffs(intvec(7*8*121,7*8));"";
+   ring r = 0,(b,c,t),dp;
+   ideal I = -13b6c3t+4b5c4t,-10b4c2t-5b4ct2;
+   primecoeffs(I);
+}  
+///////////////////////////////////////////////////////////////////////////////

Singular/LIB/presolve.lib

@@ -1 +1 @@
 //last change: 13.02.2001 (Eric Westenberger)
 ///////////////////////////////////////////////////////////////////////////////
-version="$Id: presolve.lib,v 1.18 2001-08-27 14:47:57 Singular Exp $";
+version="$Id: presolve.lib,v 1.19 2001-11-05 16:06:29 pfister Exp $";
 category="Symbolic-numerical solving";
 info="
@@ -21 +21 @@
  sortvars(id[n1,p1..]); sort vars w.r.t. complexity in id [different blocks]
  valvars(id[..]);       valuation of vars w.r.t. to their complexity in id
-           (parameters in square brackets [] are optional)
+ idealsSimplify(id);    ideal I = eliminate(Id,m) m is a product of variables
+                       which are only linearly involved in the generators of id
+ idealSplit(id,tF,fS);  a list of ideals such that their intersection 
+                        has the same radical as id
+                       ( parameters in square brackets [] are optional)
 ";
 
@@ -1166 +1170 @@
 }
 ///////////////////////////////////////////////////////////////////////////////
+proc idealSplit(ideal I,list #)
+"USAGE:  idealSplit(id,timeF,timeS);  id ideal and optioonal 
+         timeF ,timeS integers to bound the time which can be used
+         for factorization resp. standard basis computation
+RETURN:  a list of ideals such that their intersection 
+         has the same radical as id
+EXAMPLE: example idealSplit; shows an example
+"
+{
+   option(redSB);
+   int j,k,e;
+   int i=1;
+   int l=attrib(I,"isSB");
+   ideal J;
+   int timeF;
+   int timeS;
+   list re,fac,te;
+
+   if(size(#)==1)
+   {
+     if(typeof(#[1])=="ideal")
+     {
+        re=#;
+     }
+     else
+     {
+       timeF=#[1];
+     }
+   }
+   if(size(#)==2)
+   {
+     if(typeof(#[1])=="list")
+     {
+        re=#[1];
+        timeF=#[2];
+     }
+     else
+     {
+       timeF=#[1];
+       timeS=#[2];
+     }
+   }
+   if(size(#)==3){re=#[1];timeF=#[2];timeS=#[3];}
+
+   fac=timeFactorize(I[1],timeF);
+
+   while((size(fac[1])==2)&&(i<size(I)))
+   { 
+      i++; 
+      fac=timeFactorize(I[i],timeF);
+   }
+   if(size(fac[1])>2)
+   {
+      for(j=2;j<=size(fac[1]);j++)
+      {
+         I[i]=fac[1][j];
+         attrib(I,"isSB",1);
+         e=1;
+         k=0;
+         while(k<size(re))
+         {
+            k++;
+            if(size(reduce(re[k],I))==0){e=0;break;}
+            attrib(re[k],"isSB",1);
+            if(size(reduce(I,re[k]))==0){re=delete(re,k);k--;}
+         }
+         if(e)
+         {
+            if(l)
+            {
+               J=I;
+               J[i]=0;
+               J=simplify(J,2);
+               attrib(J,"isSB",1);
+               re=idealSplit(std(J,fac[1][j]),re,timeF,timeS);
+            }
+            else
+            {
+               re=idealSplit(timeStd(I,timeS),re,timeF,timeS);
+            }
+         }
+      } 
+      return(re); 
+   }
+   J=timeStd(I,timeS);
+   attrib(I,"isSB",1);
+   if(size(reduce(J,I))==0){return(re+list(I));}
+   return(re+idealSplit(J,re,timeF,timeS));
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r=32003,(b,s,t,u,v,w,x,y,z),dp;
+   ideal i=
+   bv+su,
+   bw+tu,
+   sw+tv,
+   by+sx,
+   bz+tx,
+   sz+ty,
+   uy+vx,
+   uz+wx,
+   vz+wy,
+   bvz;
+   idealSplit(i);
+}
+///////////////////////////////////////////////////////////////////////////////
+proc idealSimplify(ideal J,list #)
+"USAGE:  idealSimplify(id);  id ideal
+RETURN:  ideal I = eliminate(Id,m) m is a product of variables
+         which are only linearly involved in the generators of id
+EXAMPLE: example idealSimplify; shows an example
+"
+{
+   ideal I=J;
+   if(size(#)!=0){I=#[1];}
+   def R=basering;
+   matrix M=jacob(I);
+   ideal ma=maxideal(1);
+   int i,j,k;
+   map phi;
+
+   for(i=1;i<=nrows(M);i++)
+   {
+      for(j=1;j<=ncols(M);j++)
+      {       
+         if(deg(M[i,j])==0)
+         {
+            ma[j]=(-1/M[i,j])*(I[i]-M[i,j]*var(j));
+            phi=R,ma;
+            I=phi(I);
+            J=phi(J);
+            for(k=1;k<=ncols(I);k++){I[k]=cleardenom(I[k]);}
+            M=jacob(I);
+         }
+      }
+   }
+   J=simplify(J,2);
+   for(i=1;i<=size(J);i++){J[i]=cleardenom(J[i]);}
+   return(J);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r=0,(x,y,z,w,t),dp;
+   ideal i=
+   t,
+   x3+y2+2z,
+   x2+3y,
+   x2+y2+z2,
+   w2+z;
+   ideal j=idealSimplify(i);
+   ideal k=eliminate(i,zyt);
+   reduce(k,std(j));
+   reduce(j,std(k));
+}
+
+///////////////////////////////////////////////////////////////////////////////
+
 /*
 

Singular/LIB/standard.lib

@@ -1 +1 @@
 //////////////////////////////////////////////////////////////////////////////
-version="$Id: standard.lib,v 1.60 2001-08-28 12:54:05 Singular Exp $";
+version="$Id: standard.lib,v 1.61 2001-11-05 16:05:08 pfister Exp $";
 category="Miscellaneous";
 info="
@@ -14 +14 @@
  fprintf(link,fmt,..)   writes formatted string to link
  printf(fmt,...)        displays formatted string
+ timeStd(i,d)           std(i) if the standard basis computation finished after
+                        d-1 seconds and i otherwhise
+ timeFactorize(p,d)     factorize(p) if the factorization finished after d-1 
+                        seconds otherwhise f is considered to be irreducible
+factorH(p)              changes variables to become the last variable the 
+                        principal one in the multivariate factorization and 
+                        factorizes then the polynomial
+
 ";
 
@@ -1025 +1033 @@
 }
 
+proc timeFactorize(poly i,list #)
+"USAGE:  timeFactorize(p,d)  poly p , integer d
+RETURN:  factorize(p) if the factorization finished after d-1 
+         seconds otherwhise f is considered to be irreducible
+EXAMPLE: example timeFactorize; shows an example
+"
+{
+  def P=basering;
+  if (size(#) > 0)
+  {
+    if (system("with", "MP"))
+    {
+      if ((typeof(#[1]) == "int")&&(#[1]))
+      {
+        int wait = #[1];
+        int j = 10;
+
+        string bs = nameof(basering);
+        link l_fork = "MPtcp:fork";
+        open(l_fork);
+        write(l_fork, quote(system("pid")));
+        int pid = read(l_fork);
+        write(l_fork, quote(timeFactorize(eval(i))));
+
+        // sleep in small intervalls for appr. one second
+        if (wait > 0)
+        {
+          while(j < 1000000)
+          {
+            if (status(l_fork, "read", "ready", j)) {break;}
+            j = j + j;
+          }
+        }
+
+        // sleep in intervalls of one second from now on
+        j = 1;
+        while (j < wait)
+        {
+          if (status(l_fork, "read", "ready", 1000000)) {break;}
+          j = j + 1;
+        }
+
+        if (status(l_fork, "read", "ready"))
+        {
+          def result = read(l_fork);
+          if (bs != nameof(basering))
+          {
+            def PP = basering;
+            setring P;
+            def result = imap(PP, result);
+            kill PP;
+          }
+          kill (l_fork);
+        }
+        else
+        {
+          list result;
+          intvec v=1,1;
+          result[1]=list(1,i);
+          result[2]=v;
+          j = system("sh", "kill " + string(pid));
+        }
+        return (result);
+      }
+    }
+  }
+  return(factorH(i));
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r=0,(x,y),dp;
+   poly p=((x2+y3)^2+xy6)*((x3+y2)^2+x10y);
+   p=p^2;
+   timeFactorize(p,2);
+   timeFactorize(p,20);
+}
+
+proc timeStd(ideal i,list #)
+"USAGE:  timeStd(i,d), i ideal, d integer
+RETURN:  std(i) if the standard basis computation finished after
+         d-1 seconds and i otherwhise
+EXAMPLE: example timeStd; shows an example
+"
+{
+  def P=basering;
+  if (size(#) > 0)
+  {
+    if (system("with", "MP"))
+    {
+      if ((typeof(#[1]) == "int")&&(#[1]))
+      {
+        int wait = #[1];
+        int j = 10;
+
+        string bs = nameof(basering);
+        link l_fork = "MPtcp:fork";
+        open(l_fork);
+        write(l_fork, quote(system("pid")));
+        int pid = read(l_fork);
+        write(l_fork, quote(timeStd(eval(i))));
+
+        // sleep in small intervalls for appr. one second
+        if (wait > 0)
+        {
+          while(j < 1000000)
+          {
+            if (status(l_fork, "read", "ready", j)) {break;}
+            j = j + j;
+          }
+        }
+        j = 1;
+        while (j < wait)
+        {
+          if (status(l_fork, "read", "ready", 1000000)) {break;}
+          j = j + 1;
+        }
+        if (status(l_fork, "read", "ready"))
+        {
+          def result = read(l_fork);
+          if (bs != nameof(basering))
+          {
+            def PP = basering;
+            setring P;
+            def result = imap(PP, result);
+            kill PP;
+          }
+          kill (l_fork);
+        }
+        else
+        {
+          ideal result=i;
+          j = system("sh", "kill " + string(pid));
+        }
+        return (result);
+      }
+    }
+  }
+  return(std(i));
+}
+example
+{ "EXAMPLE:"; echo = 2;
+   ring r=32003,(a,b,c,d,e),dp;
+   int n=6;
+   ideal i=
+   a^n-b^n,
+   b^n-c^n,
+   c^n-d^n,
+   d^n-e^n,
+   a^(n-1)*b+b^(n-1)*c+c^(n-1)*d+d^(n-1)*e+e^(n-1)*a;
+   timeStd(i,2);
+   timeStd(i,20);
+}
+
+proc factorH(poly p)
+"USAGE:  factorH(p)  p poly
+RETURN:  factorize(p)
+NOTE:    changes variables to become the last variable the principal
+         one in the multivariate factorization and factorizes then
+         the polynomial
+EXAMPLE: example factorH; shows an example
+"
+{
+   def R=basering;
+   int i,j;
+   int n=1;
+   int d=nrows(coeffs(p,var(1)));
+   for(i=1;i<=nvars(R);i++)
+   {
+      j=nrows(coeffs(p,var(i)));
+      if(d>j)
+      {
+         n=i;
+         d=j;
+      }
+   }
+   ideal ma=maxideal(1);    //die letzte Variable ist die Hauptvariable
+   ma[nvars(R)]=var(n);
+   ma[n]=var(nvars(R));
+   map phi=R,ma;
+   list fac=factorize(phi(p));
+   list re=phi(fac);
+   return(re);
+}
+example
+{ "EXAMPLE:"; echo = 2;
+  system("random",992851144);
+  ring r=32003,(x,y,z,w,t),lp;
+  poly p=y2w9+yz7t-yz5w4-z2w4t4-w8t3;
+  factorize(p);  //fast
+  system("random",992851262);
+  factorize(p);  //slow
+  system("random",992851262);
+  factorH(p);
+}
+
 /*
 proc minres(list #)