Changeset f11d7b in git

Timestamp:

Jan 19, 2004, 12:26:20 PM (19 years ago)

Author:

Hans Schönemann <hannes@…>

Branches:

(u'spielwiese', '828514cf6e480e4bafc26df99217bf2a1ed1ef45')

Children:

b44d16e5c775914df07cc382ca039bb32caae77a

Parents:

bc9f64c8ededf224f4b2323f83f51d42ba00716e

Message:

*hannes: gcd via NTL


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

Location:

factory

Files:

• NTLconvert.cc (modified) (9 diffs)
• NTLconvert.h (modified) (3 diffs)
• cf_gcd.cc (modified) (3 diffs)

factory/NTLconvert.cc

@@ -1 @@
-/* $Id: NTLconvert.cc,v 1.12 2003-10-15 17:19:38 Singular Exp $ */
+/* $Id: NTLconvert.cc,v 1.13 2004-01-19 11:26:19 Singular Exp $ */
 #include <config.h>
 
@@ -20 +20 @@
 #include <NTL/ZZXFactoring.h>
 #include <NTL/ZZ_pXFactoring.h>
+#include <NTL/lzz_pXFactoring.h>
 #include <NTL/GF2XFactoring.h>
 #include <NTL/ZZ_pEXFactoring.h>
+#include <NTL/lzz_pEXFactoring.h>
 #include <NTL/GF2EXFactoring.h>
 #include <NTL/tools.h>
@@ -146 +148 @@
   return ntl_poly;
 }
+zz_pX convertFacCF2NTLzzpX(CanonicalForm f)
+{
+  zz_pX ntl_poly;
+
+  CFIterator i;
+  i=f;
+
+  int j=0;
+  int NTLcurrentExp=i.exp();
+  int largestExp=i.exp();
+  int k;
+
+  // we now build up the NTL-polynomial
+  ntl_poly.SetMaxLength(largestExp+1);
+
+  for (;i.hasTerms();i++)
+  {
+    for (k=NTLcurrentExp;k>i.exp();k--)
+    {
+      SetCoeff(ntl_poly,k,0);
+    }
+    NTLcurrentExp=i.exp();
+
+    CanonicalForm c=i.coeff();
+    if (!c.isImm()) c.mapinto(); //c%= getCharacteristic();
+    if (!c.isImm())
+    {  //This case will never happen if the characteristic is in fact a prime
+       // number, since all coefficients are represented as immediates
+       #ifndef NOSTREAMIO
+       cout<<"convertFacCF2NTLzz_pX: coefficient not immediate! : "<<f<<"\n";
+       #else
+       printf("convertFacCF2NTLzz_pX: coefficient not immediate!, char=%d\n",
+              getCharacteristic());
+       #endif
+       exit(1);
+    }
+    else
+    {
+      SetCoeff(ntl_poly,NTLcurrentExp,c.intval());
+    }
+    NTLcurrentExp--;
+  }
+
+  //Set the remaining coefficients of ntl_poly to zero.
+  // This is necessary, because NTL internally
+  // also stores powers with zero coefficient,
+  // whereas factory stores tuples of degree and coefficient
+  //leaving out tuples if the coefficient equals zero
+  for (k=NTLcurrentExp;k>=0;k--)
+  {
+    SetCoeff(ntl_poly,k,0);
+  }
+
+  //normalize the polynomial and return it
+  ntl_poly.normalize();
+
+  return ntl_poly;
+}
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -261 +321 @@
 }
 
+CanonicalForm convertNTLzzpX2CF(zz_pX poly,Variable x)
+{
+  //printf("convertNTLzzpX2CF\n");
+  CanonicalForm bigone;
+
+
+  if (deg(poly)>0)
+  {
+    // poly is non-constant
+    bigone=0;
+    bigone.mapinto();
+    // Compute the canonicalform coefficient by coefficient,
+    // bigone summarizes the result.
+    for (int j=0;j<deg(poly)+1;j++)
+    {
+      if (coeff(poly,j)!=0)
+      {
+        bigone+=(power(x,j)*CanonicalForm(to_long(rep(coeff(poly,j)))));
+      }
+    }
+  }
+  else
+  {
+    // poly is immediate
+    bigone=CanonicalForm(to_long(rep(coeff(poly,0))));
+    bigone.mapinto();
+  }
+  return bigone;
+}
+
+CanonicalForm convertNTLZZX2CF(ZZX polynom,Variable x)
+{
+  //printf("convertNTLZZX2CF\n");
+  CanonicalForm bigone;
+
+  // Go through the vector e and build up the CFFList
+  // As usual bigone summarizes the result
+  bigone=0;
+  ZZ coefficient;
+
+  for (int j=0;j<=deg(polynom);j++)
+  {
+    coefficient=coeff(polynom,j);
+    if (!IsZero(coefficient))
+    {
+      bigone += (power(x,j)*convertZZ2CF(coefficient));
+    }
+  }
+  return bigone;
+}
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -357 +467 @@
   {
     rueckgabe.append(CFFactor(convertNTLZZpX2CF(e[i].a,x),e[i].b));
+  }
+  if(isOn(SW_USE_NTL_SORT)) rueckgabe.sort(NTLcmpCF);
+  // the multiplicity at pos 1
+  if (!IsOne(multi))
+    rueckgabe.insert(CFFactor(CanonicalForm(to_long(rep(multi))),1));
+  return rueckgabe;
+}
+CFFList convertNTLvec_pair_zzpX_long2FacCFFList
+                                  (vec_pair_zz_pX_long e,zz_p multi,Variable x)
+{
+  //printf("convertNTLvec_pair_zzpX_long2FacCFFList\n");
+  CFFList rueckgabe;
+  zz_pX polynom;
+  long exponent;
+  CanonicalForm bigone;
+
+  // Maybe, e may additionally be sorted with respect to increasing degree of x
+  // but this is not
+  //important for the factorization, but nevertheless would take computing time,
+  // so it is omitted
+
+
+  // Go through the vector e and compute the CFFList
+  // again bigone summarizes the result
+  for (int i=e.length()-1;i>=0;i--)
+  {
+    rueckgabe.append(CFFactor(convertNTLzzpX2CF(e[i].a,x),e[i].b));
   }
   if(isOn(SW_USE_NTL_SORT)) rueckgabe.sort(NTLcmpCF);
@@ -629 +766 @@
   CanonicalForm bigone;
 
-  // Maybe, e may additionally be sorted with respect to increasing degree of x, but this is not
-  //important for the factorization, but nevertheless would take computing time, so it is omitted
-
-
   // Go through the vector e and build up the CFFList
   // As usual bigone summarizes the result
   for (int i=e.length()-1;i>=0;i--)
   {
-    bigone=0;
     ZZ coefficient;
     polynom=e[i].a;
     exponent=e[i].b;
-    long coeff_long;
-
-    for (int j=0;j<deg(polynom)+1;j++)
-    {
-      coefficient=coeff(polynom,j);
-      if (!IsZero(coefficient))
-      {
-        bigone += (power(x,j)*convertZZ2CF(coefficient));
-      }
-    }
-
+    bigone=convertNTLZZX2CF(polynom,x);
     //append the converted polynomial to the list
     rueckgabe.append(CFFactor(bigone,exponent));
@@ -684 +806 @@
 {
   return convertNTLZZpX2CF(rep(coefficient),x);
+}
+CanonicalForm convertNTLzzpE2CF(zz_pE coefficient,Variable x)
+{
+  return convertNTLzzpX2CF(rep(coefficient),x);
 }
 
@@ -751 +877 @@
   return rueckgabe;
 }
+CFFList convertNTLvec_pair_zzpEX_long2FacCFFList(vec_pair_zz_pEX_long e,zz_pE multi,Variable x,Variable alpha)
+{
+  CFFList rueckgabe;
+  zz_pEX polynom;
+  long exponent;
+  CanonicalForm bigone;
+
+  // Maybe, e may additionally be sorted with respect to increasing degree of x, but this is not
+  //important for the factorization, but nevertheless would take computing time, so it is omitted
+
+  // Go through the vector e and build up the CFFList
+  // As usual bigone summarizes the result during every loop
+  for (int i=e.length()-1;i>=0;i--)
+  {
+    bigone=0;
+
+    polynom=e[i].a;
+    exponent=e[i].b;
+
+    for (int j=0;j<deg(polynom)+1;j++)
+    {
+      if (IsOne(coeff(polynom,j)))
+      {
+        bigone+=power(x,j);
+      }
+      else
+      {
+        CanonicalForm coefficient=convertNTLzzpE2CF(coeff(polynom,j),alpha);
+        if (coeff(polynom,j)!=0)
+        {
+          bigone += (power(x,j)*coefficient);
+        }
+      }
+    }
+    //append the computed polynomials together with its exponent to the CFFList
+    rueckgabe.append(CFFactor(bigone,exponent));
+  }
+  if(isOn(SW_USE_NTL_SORT)) rueckgabe.sort(NTLcmpCF);
+  // Start by appending the multiplicity
+  if (!IsOne(multi))
+    rueckgabe.insert(CFFactor(convertNTLzzpE2CF(multi,alpha),1));
+
+  //return the computed CFFList
+  return rueckgabe;
+}
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -903 +1074 @@
   return result;
 }
+zz_pEX convertFacCF2NTLzz_pEX(CanonicalForm f, zz_pX mipo)
+{
+  zz_pE::init(mipo);
+  zz_pEX result;
+  CFIterator i;
+  i=f;
+
+  int j=0;
+  int NTLcurrentExp=i.exp();
+  int largestExp=i.exp();
+  int k;
+
+  result.SetMaxLength(largestExp+1);
+  for(;i.hasTerms();i++)
+  {
+    for(k=NTLcurrentExp;k>i.exp();k--) SetCoeff(result,k,0);
+    NTLcurrentExp=i.exp();
+    CanonicalForm c=i.coeff();
+    zz_pX cc=convertFacCF2NTLzzpX(c);
+    //ZZ_pE ccc;
+    //conv(ccc,cc);
+    SetCoeff(result,NTLcurrentExp,to_zz_pE(cc));
+    NTLcurrentExp--;
+  }
+  for(k=NTLcurrentExp;k>=0;k--) SetCoeff(result,k,0);
+  result.normalize();
+  return result;
+}
 #endif

factory/NTLconvert.h

@@ -1 @@
-/* $Id: NTLconvert.h,v 1.6 2003-08-28 14:15:18 Singular Exp $ */
+/* $Id: NTLconvert.h,v 1.7 2004-01-19 11:26:20 Singular Exp $ */
 #ifndef INCL_NTLCONVERT_H
 #define INCL_NTLCONVERT_H
@@ -23 +23 @@
 #include <NTL/ZZXFactoring.h>
 #include <NTL/ZZ_pXFactoring.h>
+#include <NTL/lzz_pXFactoring.h>
 #include <NTL/GF2XFactoring.h>
 #include "int_int.h"
@@ -34 +35 @@
 
 ZZ_pX convertFacCF2NTLZZpX(CanonicalForm f);
+zz_pX convertFacCF2NTLzzpX(CanonicalForm f);
 GF2X convertFacCF2NTLGF2X(CanonicalForm f);
 CanonicalForm convertNTLZZpX2CF(ZZ_pX poly,Variable x);
+CanonicalForm convertNTLzzpX2CF(zz_pX poly,Variable x);
 CanonicalForm convertNTLGF2X2CF(GF2X poly,Variable x);
+CanonicalForm convertNTLZZX2CF(ZZX polynom,Variable x);
 CFFList convertNTLvec_pair_ZZpX_long2FacCFFList(vec_pair_ZZ_pX_long e,ZZ_p multi,Variable x);
+CFFList convertNTLvec_pair_zzpX_long2FacCFFList(vec_pair_zz_pX_long e,zz_p multi,Variable x);
 
 CFFList convertNTLvec_pair_GF2X_long2FacCFFList(vec_pair_GF2X_long e,GF2 multi,Variable x);

factory/cf_gcd.cc

@@ -1 +1 @@
 /* emacs edit mode for this file is -*- C++ -*- */
-/* $Id: cf_gcd.cc,v 1.21 2001-01-30 13:41:12 Singular Exp $ */
+/* $Id: cf_gcd.cc,v 1.22 2004-01-19 11:26:20 Singular Exp $ */
 
 #include <config.h>
@@ -18 +18 @@
 #include "ftmpl_functions.h"
 
+#ifdef HAVE_NTL
+#include "NTLconvert.h"
+#endif
+
 static CanonicalForm gcd_poly( const CanonicalForm & f, const CanonicalForm& g, bool modularflag );
 
@@ -170 +174 @@
 gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
 {
-    CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
-    int p, i, n;
-
-    if ( primitive ) {
-        f = F;
-        g = G;
-        c = 1;
-    }
-    else {
-        CanonicalForm cF = content( F ), cG = content( G );
-        f = F / cF;
-        g = G / cG;
-        c = bgcd( cF, cG );
-    }
-    cg = gcd( f.lc(), g.lc() );
-    cl = ( f.lc() / cg ) * g.lc();
+#ifdef HAVE_NTL 
+  if (isOn(SW_USE_NTL_GCD))
+  {
+    if ( getCharacteristic() > 0 )
+    {
+      zz_pContext ccc(getCharacteristic());
+      ccc.restore();
+      zz_p::init(getCharacteristic());
+      zz_pX F1=convertFacCF2NTLzzpX(F);
+      zz_pX G1=convertFacCF2NTLzzpX(G);
+      zz_pX R=GCD(F1,G1);
+      return convertNTLzzpX2CF(R,F.mvar());
+    }
+    else
+    {
+      ZZX F1=convertFacCF2NTLZZX(F);
+      ZZX G1=convertFacCF2NTLZZX(G);
+      ZZX R=GCD(F1,G1);
+      return convertNTLZZX2CF(R,F.mvar());
+    }
+  }
+#endif
+  CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
+  int p, i, n;
+
+  if ( primitive )
+  {
+    f = F;
+    g = G;
+    c = 1;
+  }
+  else
+  {
+    CanonicalForm cF = content( F ), cG = content( G );
+    f = F / cF;
+    g = G / cG;
+    c = bgcd( cF, cG );
+  }
+  cg = gcd( f.lc(), g.lc() );
+  cl = ( f.lc() / cg ) * g.lc();
 //     B = 2 * cg * tmin(
 //         maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
 //         maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
 //         )+1;
-    M = tmin( maxNorm(f), maxNorm(g) );
-    BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
-    q = 0;
-    i = cf_getNumSmallPrimes() - 1;
-    while ( true ) {
-        B = BB;
-        while ( i >= 0 && q < B ) {
-            p = cf_getSmallPrime( i );
-            i--;
-            while ( i >= 0 && mod( cl, p ) == 0 ) {
-                p = cf_getSmallPrime( i );
-                i--;
-            }
-            setCharacteristic( p );
-            Dp = gcd( mapinto( f ), mapinto( g ) );
-            Dp = ( Dp / Dp.lc() ) * mapinto( cg );
-            setCharacteristic( 0 );
-            if ( Dp.degree() == 0 ) return c;
-            if ( q.isZero() ) {
-                D = mapinto( Dp );
-                q = p;
-                B = power(CanonicalForm(2),D.degree())*M+1;
-            }
-            else {
-                if ( Dp.degree() == D.degree() ) {
-                    chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
-                    q = newq;
-                    D = newD;
-                }
-                else if ( Dp.degree() < D.degree() ) {
-                    // all previous p's are bad primes
-                    q = p;
-                    D = mapinto( Dp );
-                    B = power(CanonicalForm(2),D.degree())*M+1;
-                }
-                // else p is a bad prime
-            }
+  M = tmin( maxNorm(f), maxNorm(g) );
+  BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
+  q = 0;
+  i = cf_getNumSmallPrimes() - 1;
+  while ( true )
+  {
+    B = BB;
+    while ( i >= 0 && q < B )
+    {
+      p = cf_getSmallPrime( i );
+      i--;
+      while ( i >= 0 && mod( cl, p ) == 0 )
+      {
+        p = cf_getSmallPrime( i );
+        i--;
+      }
+      setCharacteristic( p );
+      Dp = gcd( mapinto( f ), mapinto( g ) );
+      Dp = ( Dp / Dp.lc() ) * mapinto( cg );
+      setCharacteristic( 0 );
+      if ( Dp.degree() == 0 )
+        return c;
+      if ( q.isZero() )
+      {
+        D = mapinto( Dp );
+        q = p;
+        B = power(CanonicalForm(2),D.degree())*M+1;
+      }
+      else
+      {
+        if ( Dp.degree() == D.degree() )
+        {
+          chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
+          q = newq;
+          D = newD;
         }
-        if ( i >= 0 ) {
-            // now balance D mod q
-            D = pp( balance( D, q ) );
-            if ( divides( D, f ) && divides( D, g ) )
-                return D * c;
-            else
-                q = 0;
+        else if ( Dp.degree() < D.degree() )
+        {
+          // all previous p's are bad primes
+          q = p;
+          D = mapinto( Dp );
+          B = power(CanonicalForm(2),D.degree())*M+1;
         }
-        else {
-            return gcd_poly( F, G, false );
-        }
-        DEBOUTLN( cerr, "another try ..." );
-    }
+        // else p is a bad prime
+      }
+    }
+    if ( i >= 0 )
+    {
+      // now balance D mod q
+      D = pp( balance( D, q ) );
+      if ( divides( D, f ) && divides( D, g ) )
+        return D * c;
+      else
+        q = 0;
+    }
+    else
+      return gcd_poly( F, G, false );
+    DEBOUTLN( cerr, "another try ..." );
+  }
 }