Changeset dc79bd in git for libpolys/polys/ext_fields/algext.cc

Timestamp:

Oct 4, 2012, 8:16:12 PM (12 years ago)

Author:

Oleksandr Motsak <motsak@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

20c54056198f825ef62888512342886cf96c82a6

Parents:

9f6cc091f81a4d6c75b4696ae3e84bc53e419051

git-author:

Oleksandr Motsak <motsak@mathematik.uni-kl.de>2012-10-04 20:16:12+02:00

git-committer:

Oleksandr Motsak <motsak@mathematik.uni-kl.de>2012-10-05 18:15:05+02:00

Message:

Major Update for Enumerators + Fixes for Algrbraic & Transcendental extensions

General:
chg: cleanup + documentation + additional assumes

Enumerators:
chg: some polish on Enumerators
add: CRecursivePolyCoeffsEnumerator<ConverterPolicy> for recursive treatment of converted coeffs as Enumerators

Coeffs:
chg: use mpz_lcm for readability in nlClearDenominators + cleanup
add: nlClear*NoPositiveLead variants should not make LC positive
chg: all nlClear* are not static in order to be usable from alg / trans exts.
fix: fixed a bug in ndClearContent

AlgExt:
add: nCoeff_is_Q_algext for checking an alg. ext. of Q
add: naClear* for alg. ext. over Q
NOTE: coeffs are polynomials in Q[a] - one should simply consider each of them recursively as a collection of numbers...
NOTE: compute GCDs over Alg. Ext... + gcds of (int.) numbers!?
NOTE: trying to be conform with older Singular: no negative leading coeff. normalization
chg: Alg. Ext: use singclap_*gcd (instead of Frank's gcd-stuff)

p_poly:
add: p_Cleardenom_n/p_Cleardenom also clear content afterwards...
chg: major and minor changes to p_Content/p_Cleardenom_n/p_Cleardenom + cleanup
add: additionally trying to assure positive leading coeff. after p_Content/p_Cleardenom_n(/p_Cleardenom?)
NOTE: which should not be needed as n_ClearDenominators/n_ClearContent are supposed to assure that themselves!
add: more assumes to p_polys.cc
NOTE: massive usage of enumerators form p_* causes problems - only doing that for Q_a()!
NOTE: do -normalization over Q(x...)

TransExt:
add: ntClear* for trans. ext's
fix: correct ntGetDenom/ntGetNumerator (thanks to pSubstPar),
NOTE: no negative denominator out of ntGetNumerator/ntGetDenom!
add: first inefficient ntClearContent/Q  and ntClearDenominators/Q & F_p impl.
NOTE:  careful with the use of nlClear* ! (only over Q!)
add: added ntTest to transext.cc on most in/outs + use ntInit(poly)!

NOTE: trying to fix the monic-poly-gcd problem in ntClearDenominators!

File:

• libpolys/polys/ext_fields/algext.cc (modified) (7 diffs)

libpolys/polys/ext_fields/algext.cc

@@ -44 +44 @@
 #include <polys/simpleideals.h>
 
+#include <polys/PolyEnumerator.h>
+
 #ifdef HAVE_FACTORY
+#include <factory/factory.h>
 #include <polys/clapconv.h>
-#include <factory/factory.h>
+#include <polys/clapsing.h>
 #endif
 
-#include "ext_fields/algext.h"
+   
+#include <polys/ext_fields/algext.h>
 #define TRANSEXT_PRIVATES 1
-#include "ext_fields/transext.h"
+#include <polys/ext_fields/transext.h>
 
 #ifdef LDEBUG
@@ -103 +107 @@
 void     naDelete(number *a, const coeffs cf);
 void     naCoeffWrite(const coeffs cf, BOOLEAN details);
-number   naIntDiv(number a, number b, const coeffs cf);
+//number   naIntDiv(number a, number b, const coeffs cf);
 const char * naRead(const char *s, number *a, const coeffs cf);
 
@@ -652 +656 @@
   naTest(a); naTest(b);
   if ((a == NULL) && (b == NULL)) WerrorS(nDivBy0);
-  return (number)p_Gcd((poly)a, (poly)b, naRing);
+  const ring R = naRing;
+  return (number) singclap_gcd(p_Copy((poly)a, R), p_Copy((poly)b, R), R);
+//  return (number)p_Gcd((poly)a, (poly)b, naRing);
 }
 
@@ -660 +666 @@
   if (a == NULL) WerrorS(nDivBy0);
   
-  poly aFactor = NULL; poly mFactor = NULL;
-  poly theGcd = p_ExtGcd((poly)a, aFactor, naMinpoly, mFactor, naRing);
+  poly aFactor = NULL; poly mFactor = NULL; poly theGcd = NULL;
+// singclap_extgcd!
+  const BOOLEAN ret = singclap_extgcd ((poly)a, naMinpoly, theGcd, aFactor, mFactor, naRing);
+
+  assume( !ret );
+
+//  if( ret ) theGcd = p_ExtGcd((poly)a, aFactor, naMinpoly, mFactor, naRing);
   
   naTest((number)theGcd); naTest((number)aFactor); naTest((number)mFactor);
@@ -845 +856 @@
 }
 
-static void naClearContent(ICoeffsEnumerator& /*numberCollectionEnumerator*/, number& c, const coeffs cf)
+
+static void naClearContent(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
 {
   assume(cf != NULL);
   assume(getCoeffType(cf) == ID);
-  assume(nCoeff_is_Q_a(cf)); // only over Q[a]/m(a), while the default impl. is used over Zp[a]/m(a) !
-  // all coeffs are given by integers!!!
-
-  c = n_Init(1, cf);
-  assume(FALSE); // TODO: NOT YET IMPLEMENTED!!!
-
-//   numberCollectionEnumerator.Reset();
-// 
-//   c = numberCollectionEnumerator.Current();
-// 
-//   n_Normalize(c, r);
-// 
-//   if (!n_IsOne(c, r))
-//   {    
-//     numberCollectionEnumerator.Current() = n_Init(1, r); // ???
-// 
-//     number inv = n_Invers(c, r);
-// 
-//     while( numberCollectionEnumerator.MoveNext() )
-//     {
-//       number &n = numberCollectionEnumerator.Current();
-//       n_Normalize(n, r); // ?
-//       n_InpMult(n, inv, r);
-//     }
-// 
-//     n_Delete(&inv, r);
-//   }
+  assume(nCoeff_is_Q_algext(cf)); // only over (Q[a]/m(a)), while the default impl. is used over Zp[a]/m(a) !
+
+  const ring   R = cf->extRing;
+  assume(R != NULL);
+  const coeffs Q = R->cf; 
+  assume(Q != NULL); 
+  assume(nCoeff_is_Q(Q));  
+
+  numberCollectionEnumerator.Reset();
+
+  if( !numberCollectionEnumerator.MoveNext() ) // empty zero polynomial?
+  {
+    c = n_Init(1, cf);
+    return;
+  }
+
+  naTest(numberCollectionEnumerator.Current());
+
+  // part 1, find a small candidate for gcd
+  int s1; int s=2147483647; // max. int
+
+  const BOOLEAN lc_is_pos=naGreaterZero(numberCollectionEnumerator.Current(),cf);
+
+  int normalcount = 0;
+  
+  poly cand1, cand;
+  
+  do
+  {
+    number& n = numberCollectionEnumerator.Current();
+    naNormalize(n, cf); ++normalcount;
+
+    naTest(n);
+    
+    cand1 = (poly)n;
+
+    s1 = p_Deg(cand1, R); // naSize?
+    if (s>s1)
+    {
+      cand = cand1;
+      s = s1;
+    }
+  } while (numberCollectionEnumerator.MoveNext() );
+
+//  assume( nlGreaterZero(cand,cf) ); // cand may be a negative integer!
+
+  cand = p_Copy(cand, R);
+  // part 2: compute gcd(cand,all coeffs)
+
+  numberCollectionEnumerator.Reset();
+
+  int length = 0;
+  while (numberCollectionEnumerator.MoveNext() )
+  {
+    number& n = numberCollectionEnumerator.Current();
+    ++length;
+    
+    if( (--normalcount) <= 0)
+      naNormalize(n, cf);
+    
+    naTest(n);
+
+//    p_InpGcd(cand, (poly)n, R);
+    
+    cand = singclap_gcd(cand, p_Copy((poly)n, R), R);
+
+//    cand1 = p_Gcd(cand,(poly)n, R); p_Delete(&cand, R); cand = cand1;
+    
+    assume( naGreaterZero((number)cand, cf) ); // ???
+/*
+    if(p_IsConstant(cand,R))
+    {
+      c = cand;
+
+      if(!lc_is_pos)
+      {
+        // make the leading coeff positive
+        c = nlNeg(c, cf);
+        numberCollectionEnumerator.Reset();
+
+        while (numberCollectionEnumerator.MoveNext() )
+        {
+          number& nn = numberCollectionEnumerator.Current();
+          nn = nlNeg(nn, cf);
+        }
+      }
+      return;
+    }
+*/
+    
+  } 
+
+  // part3: all coeffs = all coeffs / cand
+  if (!lc_is_pos)
+    cand = p_Neg(cand, R);
+
+  c = (number)cand; naTest(c);  
+
+  poly cInverse = (poly)naInvers(c, cf);
+  assume(cInverse != NULL); // c is non-zero divisor!?
+
+  
+  numberCollectionEnumerator.Reset();
+
+
+  while (numberCollectionEnumerator.MoveNext() )
+  {
+    number& n = numberCollectionEnumerator.Current();
+
+    assume( length > 0 );
+
+    if( --length > 0 )
+    {
+      assume( cInverse != NULL );
+      n = (number) p_Mult_q(p_Copy(cInverse, R), (poly)n, R);
+    }
+    else
+    {
+      n = (number) p_Mult_q(cInverse, (poly)n, R);
+      cInverse = NULL;
+      assume(length == 0);
+    }
+    
+    definiteReduce((poly &)n, naMinpoly, cf);    
+  }
+  
+  assume(length == 0);
+  assume(cInverse == NULL); //   p_Delete(&cInverse, R);
+
+  // Quick and dirty fix for constant content clearing... !?
+  CRecursivePolyCoeffsEnumerator<NAConverter> itr(numberCollectionEnumerator); // recursively treat the numbers as polys!
+
+  number cc;
+
+  extern void nlClearContentNoPositiveLead(ICoeffsEnumerator&, number&, const coeffs);
+
+  nlClearContentNoPositiveLead(itr, cc, Q); // TODO: get rid of (-LC) normalization!? 
+
+  // over alg. ext. of Q // takes over the input number
+  c = (number) p_Mult_nn( (poly)c, cc, R); 
+//      p_Mult_q(p_NSet(cc, R), , R);
+
+  n_Delete(&cc, Q);
+
+  // TODO: the above is not enough! need GCD's of polynomial coeffs...!
+/*
+  // old and wrong part of p_Content
+    if (rField_is_Q_a(r) && !CLEARENUMERATORS) // should not be used anymore if CLEARENUMERATORS is 1
+    {
+      // we only need special handling for alg. ext.
+      if (getCoeffType(r->cf)==n_algExt)
+      {
+        number hzz = n_Init(1, r->cf->extRing->cf);
+        p=ph;
+        while (p!=NULL)
+        { // each monom: coeff in Q_a
+          poly c_n_n=(poly)pGetCoeff(p);
+          poly c_n=c_n_n;
+          while (c_n!=NULL)
+          { // each monom: coeff in Q
+            d=n_Lcm(hzz,pGetCoeff(c_n),r->cf->extRing->cf);
+            n_Delete(&hzz,r->cf->extRing->cf);
+            hzz=d;
+            pIter(c_n);
+          }
+          pIter(p);
+        }
+        // hzz contains the 1/lcm of all denominators in c_n_n
+        h=n_Invers(hzz,r->cf->extRing->cf);
+        n_Delete(&hzz,r->cf->extRing->cf);
+        n_Normalize(h,r->cf->extRing->cf);
+        if(!n_IsOne(h,r->cf->extRing->cf))
+        {
+          p=ph;
+          while (p!=NULL)
+          { // each monom: coeff in Q_a
+            poly c_n=(poly)pGetCoeff(p);
+            while (c_n!=NULL)
+            { // each monom: coeff in Q
+              d=n_Mult(h,pGetCoeff(c_n),r->cf->extRing->cf);
+              n_Normalize(d,r->cf->extRing->cf);
+              n_Delete(&pGetCoeff(c_n),r->cf->extRing->cf);
+              pGetCoeff(c_n)=d;
+              pIter(c_n);
+            }
+            pIter(p);
+          }
+        }
+        n_Delete(&h,r->cf->extRing->cf);
+      }
+    }
+*/  
+
+  
+//  c = n_Init(1, cf); assume(FALSE); // TODO: NOT YET IMPLEMENTED!!!
+}
+
+
+static void naClearDenominators(ICoeffsEnumerator& numberCollectionEnumerator, number& c, const coeffs cf)
+{
+  assume(cf != NULL);
+  assume(getCoeffType(cf) == ID);
+  assume(nCoeff_is_Q_algext(cf)); // only over (Q[a]/m(a)), while the default impl. is used over Zp[a]/m(a) !
+
+  assume(cf->extRing != NULL);
+  const coeffs Q = cf->extRing->cf; 
+  assume(Q != NULL); 
+  assume(nCoeff_is_Q(Q));  
+  number n;
+  CRecursivePolyCoeffsEnumerator<NAConverter> itr(numberCollectionEnumerator); // recursively treat the numbers as polys!
+
+  extern void nlClearDenominatorsNoPositiveLead(ICoeffsEnumerator&, number&, const coeffs);
+
+  nlClearDenominatorsNoPositiveLead(itr, n, Q); // this should probably be fine...
+  c = (number)p_NSet(n, cf->extRing); // over alg. ext. of Q // takes over the input number
 }
 
@@ -951 +1151 @@
   cf->nCoeffIsEqual  = naCoeffIsEqual;
   cf->cfInvers       = naInvers;
-  cf->cfIntDiv       = naDiv;
+  cf->cfIntDiv       = naDiv; // ???
 #ifdef HAVE_FACTORY
   cf->convFactoryNSingN=naConvFactoryNSingN;
@@ -963 +1163 @@
 
   if( nCoeff_is_Q(R->cf) )
+  {
     cf->cfClearContent = naClearContent;
+    cf->cfClearDenominators = naClearDenominators;
+  }
   
   return FALSE;