Changeset 359d742 in git

Timestamp:

Jun 24, 2008, 2:52:44 PM (16 years ago)

Author:

Hans Schönemann <hannes@…>

Branches:

(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b4f17ed1d25f93d46dbe29e4b499baecc2fd51bb')

Children:

18c7a3e40a89b378bb7b4666e1f16148b663f9e8

Parents:

863f53819a571f2381b1eac36056dff30adc1904

Message:

*hannens: new versions:EZGCD_P,QGCD


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

Location:

factory

Files:

• algext.cc (modified) (5 diffs)
• algext.h (modified) (1 diff)
• cf_defs.h (modified) (2 diffs)
• cf_switches.cc (modified) (2 diffs)
• cf_switches.h (modified) (2 diffs)

factory/algext.cc

@@ -9 +9 @@
 #include "cf_algorithm.h"
 #include "algext.h"
-
+#include "fieldGCD.h"
+#include "cf_map.h"
+#include "cf_generator.h"
 
 void tryEuclid( const CanonicalForm & A, const CanonicalForm & B, const CanonicalForm M, CanonicalForm & result, bool & fail )
@@ -61 +63 @@
   Variable x = Variable(1);
   if(!extgcd( replacevar( F, a, x ), replacevar( M, a, x ), inv, b ).isOne())
+  {
+    printf("failed ");
     fail = true;
+  }
   else
     inv = replacevar( inv, a, x); // change back to alg var
@@ -84 +89 @@
 }
 
-
-CanonicalForm QGCD( const CanonicalForm & F, const CanonicalForm & G )
-{
+CanonicalForm univarQGCD( const CanonicalForm & F, const CanonicalForm & G )
+{ // F,G in Q(a)[x]
   CanonicalForm f, g, tmp, M, q, D, Dp, cl, newD, newq;
   int p, dp_deg, bound, i;
@@ -133 +137 @@
     dp_deg = degree(Dp);
 
-    if( !dp_deg ) // early termination
+    if( dp_deg == 0 ) // early termination
+    {
+      CanonicalForm inv;
+      tryInvert(Dp, M, inv, fail);
+      if(fail)
+        continue;
       return CanonicalForm(1);
+    }
 
     if( dp_deg > bound ) // current prime unlucky
@@ -171 +181 @@
   return D;
 }
+
+
+CanonicalForm QGCD( const CanonicalForm & F, const CanonicalForm & G )
+{ // F,G in Q(a)[x1,...,xn]
+  if(F.isZero())
+  {
+    if(G.isZero())
+      return G; // G is zero
+    if(G.inCoeffDomain())
+      return CanonicalForm(1);
+    return G/Lc(G); // return monic G
+  }
+  if(G.isZero()) // F is non-zero
+  {
+    if(F.inCoeffDomain())
+      return CanonicalForm(1);
+    return F/Lc(F); // return monic F
+  }
+  if(F.inCoeffDomain() || G.inCoeffDomain())
+    return CanonicalForm(1);
+
+  CanonicalForm D;
+  if (getCharacteristic()==0)
+  {
+    CanonicalForm f,g,tmp, M, q, Dp, cl, newD, newq;
+    int p, i;
+    int *bound, *other; // degree vectors
+    bool fail;
+    On(SW_RATIONAL);
+    f = F * bCommonDen(F);
+    g = G * bCommonDen(G);
+    Variable a,b;
+    if( !hasFirstAlgVar( f, a ) && !hasFirstAlgVar( g, b ))
+    {
+      // F and G are in Q[x1,...,xn], call...
+      return gcd_poly_0( f, g );
+    }
+    if( b.level() > a.level() )
+      a = b;
+    // here: a is the biggest alg. var in f and g
+    tmp = getMipo(a);
+    M = tmp * bCommonDen(tmp);
+    Off( SW_RATIONAL );
+    // here: f, g in Z[y][x1,...,xn], M in Z[y] not necessarily monic
+    // calculate upper bound for degree of gcd
+    int mv = f.level();
+    if(g.level() > mv)
+      mv = g.level();
+    // here: mv is level of the largest variable in f, g
+    bound = new int[mv+1];
+    other = new int[mv+1];
+    for(int i=1; i<=mv; i++)
+      bound[i] = other[i] = 0;
+    bound = leadDeg(f,bound);
+    other = leadDeg(g,other);
+    for(int i=1; i<=mv; i++) // entry at i=0 is not visited
+      if(other[i] < bound[i])
+        bound[i] = other[i];
+    // now bound points on the smaller vector
+    cl = lc(M) * lc(f) * lc(g);
+    q = 1;
+    D = 0;
+    for(int i=cf_getNumBigPrimes()-1; i>-1; i--)
+    {
+      p = cf_getBigPrime(i);
+      if( mod( cl, p ).isZero() ) // skip lc-bad primes
+        continue;
+
+      fail = false;
+      setCharacteristic(p);
+      tryBrownGCD( mapinto(f), mapinto(g), mapinto(M), Dp, fail );
+      setCharacteristic(0);
+      if( fail ) // M splits in char p
+        continue;
+
+      for(int i=1; i<=mv; i++)
+        other[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
+      other = leadDeg(Dp,other);
+
+      if( other==0 ) // early termination
+      {
+        // Dp is in coeff domain
+        CanonicalForm inv;
+        tryInvert(Dp,M,inv,fail); // check if zero-divisor
+        if(fail)
+          continue;
+        return CanonicalForm(1);
+      }
+
+      if( isLess(bound, other, 1, mv) ) // current prime unlucky
+        continue;
+
+      if( isLess(other, bound, 1, mv) ) // all previous primes unlucky
+      {
+        q = p;
+        D = mapinto(Dp); // shortcut CRA
+        for(int i=1; i<=mv; i++) // tighten bound
+          bound[i] = other[i];
+        continue;
+      }
+      chineseRemainder( D, q, mapinto(Dp), p, newD, newq );
+      // newD = Dp mod p
+      // newD = D mod q
+      // newq = p*q
+      q = newq;
+      if( D != newD )
+      {
+        D = newD;
+        continue;
+      }
+      On( SW_RATIONAL );
+      tmp = Farey( D, q );
+      if( fdivides( tmp, f ) && fdivides( tmp, g )) // trial division
+      {
+        Off( SW_RATIONAL );
+        return tmp;
+      }
+      Off( SW_RATIONAL );
+    }
+  }
+  // hopefully, we never reach this point
+  Off( SW_USE_QGCD );
+  D = gcd( F, G );
+  On( SW_USE_QGCD );
+  return D;
+}
+
+void tryBrownGCD( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, bool & fail )
+{// assume F,G are multivariate polys over Z/p(a) for big prime p, M "univariate" polynomial in an algebraic variable
+  printf("Brown ");
+  if(F.isZero())
+  {
+    if(G.isZero())
+    {
+      result = G; // G is zero
+      return;
+    }
+    if(G.inCoeffDomain())
+    {
+      tryInvert(G,M,result,fail);
+      return;
+    }
+    // try to make G monic modulo M
+    CanonicalForm inv;
+    tryInvert(Lc(G),M,inv,fail);
+    if(fail)
+      return;
+    result = inv*G;
+    return;
+  }
+  if(G.isZero()) // F is non-zero
+  {
+    if(F.inCoeffDomain())
+    {
+      tryInvert(F,M,result,fail);
+      return;
+    }
+    // try to make F monic modulo M
+    CanonicalForm inv;
+    tryInvert(Lc(F),M,inv,fail);
+    if(fail)
+      return;
+    result = inv*F;
+    return;
+  }
+  if(F.inCoeffDomain())
+  {
+    tryInvert(F,M,result,fail);
+    return;
+  }
+  if(G.inCoeffDomain())
+  {
+    tryInvert(G,M,result,fail);
+    return;
+  }
+  CFMap MM,NN;
+  CFArray ps(1,2);
+  ps[1] = F;
+  ps[2] = G;
+  compress(ps,MM,NN); // maps MM, NN are created
+  CanonicalForm f=MM(F);
+  CanonicalForm g=MM(G);
+  // here: f, g are compressed
+  // compute largest variable in f or g (least one is Variable(1))
+  int mv = f.level();
+  if(g.level() > mv)
+    mv = g.level();
+  // here: mv is level of the largest variable in f, g
+  if(mv == 1) // f,g univariate
+  {
+    tryEuclid(f,g,M,result,fail);
+    if(fail)
+      return;
+    result = NN(result); // do not forget to map back
+    return;
+  }
+  // here: mv > 1
+  CanonicalForm cf = vcontent(f, Variable(2)); // cf is univariate poly in var(1)
+  CanonicalForm cg = vcontent(g, Variable(2));
+  CanonicalForm c;
+  tryEuclid(cf,cg,M,c,fail);
+  if(fail)
+    return;
+  f/=cf;
+  g/=cg;
+  if(f.inCoeffDomain())
+  {
+    tryInvert(f,M,result,fail);
+    if(fail)
+      return;
+    result = NN(result);
+    return;
+  }
+  if(g.inCoeffDomain())
+  {
+    tryInvert(g,M,result,fail);
+    if(fail)
+      return;
+    result = NN(result);
+    return;
+  }
+  int *L = new int[mv+1]; // L is addressed by i from 2 to mv
+  int *N = new int[mv+1];
+  for(int i=2; i<=mv; i++)
+    L[i] = N[i] = 0;
+  L = leadDeg(f, L);
+  N = leadDeg(g, N);
+  CanonicalForm gamma;
+  tryEuclid( firstLC(f), firstLC(g), M, gamma, fail );
+  if(fail)
+    return;
+  for(int i=2; i<=mv; i++) // entry at i=1 is not visited
+    if(N[i] < L[i])
+      L[i] = N[i];
+  // L is now upper bound for degrees of gcd
+  int *dg_im = new int[mv+1]; // for the degree vector of the image we don't need any entry at i=1
+  for(int i=2; i<=mv; i++)
+    dg_im[i] = 0; // initialize
+  CanonicalForm gamma_image, m=1;
+  CanonicalForm gm=0;
+  CanonicalForm g_image, alpha, gnew, mnew;
+  FFGenerator gen = FFGenerator();
+  for(FFGenerator gen = FFGenerator(); gen.hasItems(); gen.next())
+  {
+    alpha = gen.item();
+    gamma_image = gamma(alpha, Variable(1)); // plug in alpha for var(1)
+    if(gamma_image.isZero()) // skip lc-bad points var(1)-alpha
+      continue;
+    tryBrownGCD( f(alpha, Variable(1)), g(alpha, Variable(1)), M, g_image, fail ); // recursive call with one var less
+    if(fail)
+      return;
+    if(g_image.inCoeffDomain()) // early termination
+    {
+      tryInvert(g_image,M,result,fail);
+      if(fail)
+        return;
+      result = NN(c);
+      return;
+    }
+    for(int i=2; i<=mv; i++)
+      dg_im[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
+    dg_im = leadDeg(g_image, dg_im); // dg_im cannot be NIL-pointer
+    if(isEqual(dg_im, L, 2, mv))
+    {
+      g_image /= lc(g_image); // make g_image monic over Z/p
+      g_image *= gamma_image; // multiply by multiple of image lc(gcd)
+      tryCRA( g_image, Variable(1)-alpha, gm, m, gnew, mnew, fail );
+      // gnew = gm mod m
+      // gnew = g_image mod var(1)-alpha
+      // mnew = m * (var(1)-alpha)
+      if(fail)
+        return;
+      m = mnew;
+      if(gnew == gm) // gnew did not change
+      {
+        g_image = gm / vcontent(gm, Variable(2));
+        if(fdivides(g_image,f) && fdivides(g_image,g)) // trial division
+        {
+          result = NN(c*g_image);
+          return;
+        }
+      }
+      gm = gnew;
+      continue;
+    }
+
+    if(isLess(L, dg_im, 2, mv)) // dg_im > L --> current point unlucky
+      continue;
+
+    if(isLess(dg_im, L, 2, mv)) // dg_im < L --> all previous points were unlucky
+    {
+      m = CanonicalForm(1); // reset
+      gm = 0; // reset
+      for(int i=2; i<=mv; i++) // tighten bound
+        L[i] = dg_im[i];
+    }
+  }
+  // we are out of evaluation points
+  fail = true;
+}
+
+
+void tryCRA( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew, bool & fail )
+{ // polys of level <= 1 are considered coefficients. q1,q2 are assumed to be coprime
+  // xnew = x1 mod q1 (coefficientwise in the above sense)
+  // xnew = x2 mod q2
+  // qnew = q1*q2
+  CanonicalForm tmp;
+  if(x1.level() <= 1 && x2.level() <= 1) // base case
+  {
+    tryExtgcd(q1,q2,tmp,xnew,qnew,fail);
+    if(fail)
+      return;
+    xnew = x1 + (x2-x1) * xnew * q1;
+    qnew = q1*q2;
+    xnew = mod(xnew,qnew);
+    return;
+  }
+  CanonicalForm tmp2;
+  xnew = 0;
+  qnew = q1 * q2;
+  // here: x1.level() > 1 || x2.level() > 1
+  if(x1.level() > x2.level())
+  {
+    for(CFIterator i=x1; i.hasTerms(); i++)
+    {
+      if(i.exp() == 0) // const. term
+      {
+        tryCRA(i.coeff(),q1,x2,q2,tmp,tmp2,fail);
+        if(fail)
+          return;
+        xnew += tmp;
+      }
+      else
+      {
+        tryCRA(i.coeff(),q1,0,q2,tmp,tmp2,fail);
+        if(fail)
+          return;
+        xnew += tmp * power(x1.mvar(),i.exp());
+      }
+    }
+    return;
+  }
+  // here: x1.level() <= x2.level() && ( x1.level() > 1 || x2.level() > 1 )
+  if(x2.level() > x1.level())
+  {
+    for(CFIterator j=x2; j.hasTerms(); j++)
+    {
+      if(j.exp() == 0) // const. term
+      {
+        tryCRA(x1,q1,j.coeff(),q2,tmp,tmp2,fail);
+        if(fail)
+          return;
+        xnew += tmp;
+      }
+      else
+      {
+        tryCRA(0,q1,j.coeff(),q2,tmp,tmp2,fail);
+        if(fail)
+          return;
+        xnew += tmp * power(x2.mvar(),j.exp());
+      }
+    }
+    return;
+  }
+  // here: x1.level() == x2.level() && x1.level() > 1 && x2.level() > 1
+  CFIterator i = x1;
+  CFIterator j = x2;
+  while(i.hasTerms() || j.hasTerms())
+  {
+    if(i.hasTerms())
+    {
+      if(j.hasTerms())
+      {
+        if(i.exp() == j.exp())
+        {
+          tryCRA(i.coeff(),q1,j.coeff(),q2,tmp,tmp2,fail);
+          if(fail)
+            return;
+          xnew += tmp * power(x1.mvar(),i.exp());
+          i++; j++;
+        }
+        else
+        {
+          if(i.exp() < j.exp())
+          {
+            tryCRA(i.coeff(),q1,0,q2,tmp,tmp2,fail);
+            if(fail)
+              return;
+            xnew += tmp * power(x1.mvar(),i.exp());
+            i++;
+          }
+          else // i.exp() > j.exp()
+          {
+            tryCRA(0,q1,j.coeff(),q2,tmp,tmp2,fail);
+            if(fail)
+              return;
+            xnew += tmp * power(x1.mvar(),j.exp());
+            j++;
+          }
+        }
+      }
+      else // j is out of terms
+      {
+        tryCRA(i.coeff(),q1,0,q2,tmp,tmp2,fail);
+        if(fail)
+          return;
+        xnew += tmp * power(x1.mvar(),i.exp());
+        i++;
+      }
+    }
+    else // i is out of terms
+    {
+      tryCRA(0,q1,j.coeff(),q2,tmp,tmp2,fail);
+      if(fail)
+        return;
+      xnew += tmp * power(x1.mvar(),j.exp());
+      j++;
+    }
+  }
+}
+
+
+void tryExtgcd( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & result, CanonicalForm & s, CanonicalForm & t, bool & fail )
+{
+  // F, G are univariate polynomials (i.e. they have exactly one polynomial variable)
+  // F and G must have the same level AND level > 0
+  // we try to calculate gcd(f,g) = s*F + t*G
+  // if a zero divisor is encontered, 'fail' is set to one
+  Variable a, b;
+  if( !hasFirstAlgVar(F,a) && !hasFirstAlgVar(G,b) ) // note lazy evaluation
+  {
+    result = extgcd( F, G, s, t ); // no zero divisors possible
+    return;
+  }
+  if( b.level() > a.level() )
+    a = b;
+  // here: a is the biggest alg. var in F and G
+  CanonicalForm M = getMipo(a);
+  CanonicalForm P;
+  if( degree(F) > degree(G) )
+  {
+    P=F; result=G; s=0; t=1;
+  }
+  else
+  {
+    P=G; result=F; s=1; t=0;
+  }
+  CanonicalForm inv, rem, q, u, v;
+  // here: degree(P) >= degree(result)
+  while(true)
+  {
+    tryInvert( Lc(result), M, inv, fail );
+    if(fail)
+      return;
+    // here: Lc(result) is invertible modulo M
+    q = Lc(P)*inv * power( P.mvar(), degree(P)-degree(result) );
+    rem = P - q*result;
+    // here: s*F + t*G = result
+    if( rem.isZero() )
+    {
+      s*=inv;
+      t*=inv;
+      result *= inv; // monify result
+      return;
+    }
+    P=result;
+    result=rem;
+    rem=u-q*s;
+    u=s;
+    s=rem;
+    rem=v-q*t;
+    v=t;
+    t=rem;
+  }
+}

factory/algext.h

@@ -5 +5 @@
 
 CanonicalForm QGCD( const CanonicalForm &, const CanonicalForm & );
+CanonicalForm univarQGCD( const CanonicalForm & F, const CanonicalForm & G );
 void tryEuclid( const CanonicalForm &, const CanonicalForm &, const CanonicalForm, CanonicalForm &, bool & );
 void tryInvert( const CanonicalForm &, const CanonicalForm &, CanonicalForm &, bool & );
 bool hasFirstAlgVar( const CanonicalForm &, Variable & );
+void tryBrownGCD( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, bool & fail );
+void tryCRA( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew, bool & fail );
+void tryExtgcd( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & result, CanonicalForm & s, CanonicalForm & t, bool & fail );
+
+
+

factory/cf_defs.h

@@ -1 +1 @@
 /* emacs edit mode for this file is -*- C++ -*- */
-/* $Id: cf_defs.h,v 1.14 2008-06-16 12:55:08 Singular Exp $ */
+/* $Id: cf_defs.h,v 1.15 2008-06-24 12:52:44 Singular Exp $ */
 
 #ifndef INCL_CF_DEFS_H
@@ -43 +43 @@
 const int SW_USE_GCD_P=14;
 const int SW_USE_QGCD=15;
+const int SW_USE_fieldGCD=16;
 //}}}
 

factory/cf_switches.cc

@@ -1 +1 @@
 /* emacs edit mode for this file is -*- C++ -*- */
-/* $Id: cf_switches.cc,v 1.10 2008-06-16 12:54:25 Singular Exp $ */
+/* $Id: cf_switches.cc,v 1.11 2008-06-24 12:52:44 Singular Exp $ */
 
 //{{{ docu
@@ -37 +37 @@
   On(SW_USE_EZGCD);
   //On(SW_USE_EZGCD_P); // still testing
-  //On(SW_USE_QGCD);
+  On(SW_USE_QGCD);
+  On(SW_USE_fieldGCD);
 }
 //}}}

factory/cf_switches.h

@@ -1 +1 @@
 /* emacs edit mode for this file is -*- C++ -*- */
-/* $Id: cf_switches.h,v 1.13 2008-06-16 12:54:05 Singular Exp $ */
+/* $Id: cf_switches.h,v 1.14 2008-06-24 12:52:44 Singular Exp $ */
 
 #ifndef INCL_CF_SWITCHES_H
@@ -19 +19 @@
 //
 //}}}
-const int CFSwitchesMax = 16;
+const int CFSwitchesMax = 17;
 //}}}