Changeset c28ecf in git

Timestamp:

May 19, 2011, 2:52:08 PM (12 years ago)

Author:

Frank Seelisch <seelisch@…>

Branches:

(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', '00e2e9c41af3fde1273eb3633f4c0c7c3db2579d')

Children:

4a2260efc22c19ac8eaf58d529af8decbb9420bb

Parents:

f0b01f1bf67efed9ec49d9df0d55f75284f5d968

git-author:

Frank Seelisch <seelisch@mathematik.uni-kl.de>2011-05-19 14:52:08+02:00

git-committer:

Mohamed Barakat <mohamed.barakat@rwth-aachen.de>2011-11-09 12:36:12+01:00

Message:

alg ext implementation now passing all coeffs and polys tests

Files:

• .gitignore (modified) (1 diff)
• libpolys/coeffs/longrat.cc (modified) (2 diffs)
• libpolys/polys/ext_fields/algext.cc (modified) (6 diffs)
• libpolys/polys/ext_fields/algext.h (modified) (1 diff)
• libpolys/polys/monomials/p_polys.cc (modified) (5 diffs)
• libpolys/polys/monomials/p_polys.h (modified) (1 diff)
• libpolys/tests/polys_test.h (modified) (4 diffs)

.gitignore

@@ -102 +102 @@
 ntl/src/gen_gmp_aux
 ntl/src/gen_lip_gmp_aux
+libpolys/tests/gftables
+libpolys/tests/p_Procs_*.so
 aclocal.m4
 autom4te.cache/

libpolys/coeffs/longrat.cc

@@ -223 +223 @@
     return FALSE;
   }
-  // only, if gmp memeory routine are from omalloc:
+  /* TODO: If next line is active, then computations in algebraic field
+           extensions over Q will throw a lot of assume violations although
+           everything is computed correctly and no seg fault appears.
+           Maybe the test is not appropriate in this case. */
   //omCheckAddrSize(a->z[0]._mp_d,a->z[0]._mp_alloc*BYTES_PER_MP_LIMB);
   if (a->z[0]._mp_alloc==0)
@@ -230 +233 @@
   if (a->s<2)
   {
-    // only, if gmp memeory routine are from omalloc:
+    /* TODO: If next line is active, then computations in algebraic field
+             extensions over Q will throw a lot of assume violations although
+             everything is computed correctly and no seg fault appears.
+             Maybe the test is not appropriate in this case. */
     //omCheckIf(omCheckAddrSize(a->n[0]._mp_d,a->n[0]._mp_alloc*BYTES_PER_MP_LIMB), return FALSE);
     if (a->z[0]._mp_alloc==0)

libpolys/polys/ext_fields/algext.cc

@@ -75 +75 @@
   assume(getCoeffType(cf) == naID);
   if (a == NULL) return TRUE;
-  p_Test((poly)a, naRing); // cannot use omCheckAddrSize(a, sizeof(*a)) here
+  p_Test((poly)a, naRing);
   assume(p_Deg((poly)a, naRing) <= p_Deg(naMinpoly, naRing));
   return TRUE;
@@ -90 +90 @@
 }
 
-void naDelete(number *a, const coeffs cf)
-{
-  if (a == NULL) return;
-  p_Delete((poly *)a, naRing);
+void naDelete(number * a, const coeffs cf)
+{
+  if (*a == NULL) return;
+  poly aAsPoly = (poly)(*a);
+  p_Delete(&aAsPoly, naRing);
   *a = NULL;
 }
@@ -339 +340 @@
 {
   #ifdef LDEBUG
-  omCheckAddr(p); omCheckAddr(reducer);
+  p_Test((poly)p, naRing);
+  p_Test((poly)reducer, naRing);
   #endif
   if (p_Deg(p, naRing) > 10 * p_Deg(reducer, naRing))
@@ -427 +429 @@
 {
   #ifdef LDEBUG
-  omCheckAddr(p); omCheckAddr(reducer);
+  p_Test((poly)p, naRing);
+  p_Test((poly)reducer, naRing);
   #endif
   p_PolyDiv(p, reducer, FALSE, naRing);
@@ -443 +446 @@
   naTest(a);
   if (a == NULL) WerrorS(nDivBy0);
-  poly aFactor; poly mFactor;
+  poly aFactor = NULL; poly mFactor = NULL;
   poly theGcd = p_ExtGcd((poly)a, aFactor, naMinpoly, mFactor, naRing);
-  /* the gcd must be one since naMinpoly is irreducible and a != NULL: */
+  /* the gcd must be 1 since naMinpoly is irreducible and a != NULL: */
   assume(naIsOne((number)theGcd, cf));      
   p_Delete(&theGcd, naRing);
   p_Delete(&mFactor, naRing);
+  /* printf("naInvers\n");
+     p_Write((poly)a, naRing);
+     p_Write(aFactor, naRing); */
   return (number)(aFactor);
 }
@@ -566 +572 @@
   
   #ifdef LDEBUG
-  omCheckAddr(naMinpoly);
+  p_Test((poly)naMinpoly, naRing);
   #endif
   

libpolys/polys/ext_fields/algext.h

@@ -69 +69 @@
 number   naLcm(number a, number b, const coeffs cf);
 number   naSize(number a, const coeffs cf);
-void     naDelete(number *a, const coeffs cf);
+void     naDelete(number * a, const coeffs cf);
 void     naCoeffWrite(const coeffs cf);
 number   naIntDiv(number a, number b, const coeffs cf);

libpolys/polys/monomials/p_polys.cc

@@ -1453 +1453 @@
   int divisorLE = p_GetExp(divisor, 1, r);
 //p_Write(p, r); p_Write(divisor, r);
-  while ((p != NULL) && (p_Deg(p, r) > p_Deg(divisor, r)))
+  while ((p != NULL) && (p_Deg(p, r) >= p_Deg(divisor, r)))
   {
     /* determine t = LT(p) / LT(divisor) */
@@ -1462 +1462 @@
     p_SetExp(t, 1, e, r);
     p_Setm(t, r);
-//printf("t\n"); p_Write(t, r);
     if (needResult) result = p_Add_q(result, p_Copy(t, r), r);
-//t = p_Mult_q(t, p_Copy(divisor, r), r); p_Write(t, r);
-//t = p_Neg(t, r); p_Write(t, r);
-//printf("p\n"); p_Write(p, r);
-//t = p_Add_q(p, t, r); p_Write(t, r);
-//p = t;
     p = p_Add_q(p, p_Neg(p_Mult_q(t, p_Copy(divisor, r), r), r), r);
-//printf("EXECUTION STOPPED\n"); while (1) { };
   }
   n_Delete(&divisorLC, r->cf);
+//printf("p_PolyDiv result:\n");
+//p_Write(result, r);
   return result;
+}
+
+/* returns NULL if p == NULL, otherwise makes p monic by dividing
+   by its leading coefficient (only done if this is not already 1);
+   this assumes that we are over a ground field so that division
+   is well-defined;
+   modifies p */
+void p_Monic(poly &p, ring r)
+{
+  if (p == NULL) return;
+  poly pp = p;
+  number lc = p_GetCoeff(p, r);
+  if (n_IsOne(lc, r->cf)) return;
+  number n = n_Init(1, r->cf);
+  p_SetCoeff(p, n, r);
+  p = pIter(p);
+  while (p != NULL)
+  {
+    number c = p_GetCoeff(p, r);
+    number n = n_Div(c, lc, r->cf);
+    n_Delete(&c, r->cf);
+    p_SetCoeff(p, n, r);
+    p = pIter(p);
+  }
+  n_Delete(&lc, r->cf);
+  p = pp;
 }
 
@@ -1481 +1502 @@
 poly p_GcdHelper(poly &p, poly &q, ring r)
 {
-  if (q == NULL) return p;
+  if (q == NULL)
+  {
+    /* We have to make p monic before we return it, so that if the
+       gcd is a unit in the ground field, we will actually return 1. */
+    p_Monic(p, r);
+    return p;
+  }
   else
   {
@@ -1511 +1538 @@
                     ring r)
 {
-//printf("p_ExtGcdHelper:\n");
   if (q == NULL)
   {
-    qFactor = NULL; pFactor = p_ISet(1, r); return p;
+    qFactor = NULL;
+    pFactor = p_ISet(1, r);
+    number n = p_GetCoeff(pFactor, r);
+    p_SetCoeff(pFactor, n_Invers(p_GetCoeff(p, r), r->cf), r);
+    n_Delete(&n, r->cf);
+    p_Monic(p, r);
+    return p;
   }
   else
   {
+//printf("p_ExtGcdHelper1:\n");
 //p_Write(p, r); p_Write(q, r);
     poly pDivQ = p_PolyDiv(p, q, TRUE, r);
-    poly ppFactor; poly qqFactor;
-    poly theGcd = p_ExtGcdHelper(q, ppFactor, p, qqFactor, r);
-    pFactor = qqFactor;
-    qFactor = p_Add_q(ppFactor,
-                      p_Neg(p_Mult_q(pDivQ, p_Copy(qqFactor, r), r), r),
+    poly ppFactor = NULL; poly qqFactor = NULL;
+//printf("p_ExtGcdHelper2:\n");
+//p_Write(p, r); p_Write(ppFactor, r);
+//p_Write(q, r); p_Write(qqFactor, r);
+    poly theGcd = p_ExtGcdHelper(q, qqFactor, p, ppFactor, r);
+//printf("p_ExtGcdHelper3:\n");
+//p_Write(q, r); p_Write(qqFactor, r);
+//p_Write(p, r); p_Write(ppFactor, r);
+//p_Write(theGcd, r);
+    pFactor = ppFactor;
+    qFactor = p_Add_q(qqFactor,
+                      p_Neg(p_Mult_q(pDivQ, p_Copy(ppFactor, r), r), r),
                       r);
+//printf("p_ExtGcdHelper4:\n");
+//p_Write(pFactor, r); p_Write(qFactor, r);
     return theGcd;
   }
@@ -1540 +1582 @@
 poly p_ExtGcd(poly p, poly &pFactor, poly q, poly &qFactor, ring r)
 {
-//printf("p_ExtGcd:\n");
+//printf("p_ExtGcd1:\n");
 //p_Write(p, r); p_Write(q, r);
   assume((p != NULL) || (q != NULL));  
-  poly a = p; poly b = q;
-  poly aFactor = pFactor; poly bFactor = qFactor; 
+  poly a = p; poly b = q; BOOLEAN aCorrespondsToP = TRUE;
   if (p_Deg(a, r) < p_Deg(b, r))
-  { a = q; b = p; aFactor = qFactor, bFactor = pFactor; }
+  { a = q; b = p; aCorrespondsToP = FALSE; }
   a = p_Copy(a, r); b = p_Copy(b, r);
-  return p_ExtGcdHelper(a, aFactor, b, bFactor, r);
+  poly aFactor = NULL; poly bFactor = NULL;
+  poly theGcd = p_ExtGcdHelper(a, aFactor, b, bFactor, r);
+  if (aCorrespondsToP) { pFactor = aFactor; qFactor = bFactor; }
+  else                 { pFactor = bFactor; qFactor = aFactor; }
+  return theGcd;
 }
 

libpolys/polys/monomials/p_polys.h

@@ -1839 +1839 @@
 poly      p_PolyDiv(poly &p, poly divisor, BOOLEAN needResult, ring r);
 
+/* returns NULL if p == NULL, otherwise makes p monic by dividing
+   by its leading coefficient (only done if this is not already 1);
+   this assumes that we are over a ground field so that division
+   is well-defined;
+   modifies p */
+void      p_Monic(poly &p, ring r);
+
 /* assumes that p and q are univariate polynomials in r,
    mentioning the same variable;

libpolys/tests/polys_test.h

@@ -208 +208 @@
 class PolysTestSuite : public CxxTest::TestSuite 
 {
+private:
+  void TestArithCf(const coeffs r)
+  {
+    clog << ("TEST: Simple Arithmetics: ");
+    clog << endl;
+
+    number two = n_Init(2, r);
+
+    number t = n_Init(1, r);  
+    ndInpAdd(t, t, r);  
+    TS_ASSERT( n_Equal(two, t, r) );
+    n_Delete(&t, r);
+  
+    if( getCoeffType(r) == n_Q )
+    {
+      number t = n_Init(1, r);  
+      nlInpAdd(t, t, r);
+      TS_ASSERT( n_Equal(two, t, r) );
+      n_Delete(&t, r);
+    }
+  
+    const int N = 66666;
+
+    number a = n_Init(N, r);
+   
+    clog<< "a: "; PrintSized(a, r);
+
+    clog<< "two: "; PrintSized(two, r);
+
+    number aa0 = n_Init(N*2, r);
+
+    number aa = n_Add(a, a, r);
+
+    clog<< "aa = a + a: "; PrintSized(aa, r);
+  
+    number aa2 = n_Mult(a, two, r);
+
+    clog<< "aa2 = a * 2: "; PrintSized(aa2, r);
+
+    number aa1 = n_Mult(two, a, r);
+ 
+    clog<< "aa1 = 2 * a: "; PrintSized(aa1, r);
+
+    n_Delete(&a, r);
+    n_Delete(&two, r);
+
+    a = n_Sub( aa, aa1, r );
+  
+    clog<< "a = aa - aa1: "; PrintSized(a, r);
+
+    TS_ASSERT( n_IsZero(a, r) );
+
+    n_Delete(&a, r);
+
+    a = n_Sub( aa, aa2, r );
+
+    clog<< "a = aa - aa2: "; PrintSized(a, r);
+
+    TS_ASSERT( n_IsZero(a, r) );
+
+    n_Delete(&a, r);
+
+    a = n_Sub( aa1, aa2, r );
+
+    clog<< "a = aa1 - aa2: "; PrintSized(a, r);
+
+    TS_ASSERT( n_IsZero(a, r) );
+
+    n_Delete(&a, r);
+
+    TS_ASSERT( n_Equal(aa, aa1, r) );
+    TS_ASSERT( n_Equal(aa, aa2, r) );
+    TS_ASSERT( n_Equal(aa1, aa2, r) );
+  
+    TS_ASSERT( n_Equal(aa0, aa, r) );
+    TS_ASSERT( n_Equal(aa0, aa1, r) );
+    TS_ASSERT( n_Equal(aa0, aa2, r) );
+
+    n_Delete(&aa, r);
+    n_Delete(&aa1, r);
+    n_Delete(&aa2, r);
+
+    n_Delete(&aa0, r);
+
+    clog << ( " >>> TEST DONE!" );
+    clog << endl;
+  }
+  void TestSumCf(const coeffs r, const unsigned long N)
+  {
+    clog << ( _2S("TEST: sum[0..") + _2S(N) + "]: ");
+    clog << endl;
+
+    assume( N > 0 ); // just for now...
+
+    const unsigned long ssss = (N * (N+1)) / 2;
+    
+    number sum1 = n_Init(ssss, r);
+    clog<< "N*(N+1)/2 (int: " << ssss << "): "; PrintSized(sum1, r);
+
+    number s, ss, i, res;
+
+    s = n_Init(N  , r);
+    i = n_Init(N+1, r);
+    ndInpMult(s, i, r);
+    n_Delete(&i, r);
+    
+    clog<< "N*(N+1): ("<< N*(N+1) << ")"; PrintSized(s, r);  
+    
+    i = n_Init(2, r);
+    clog<< "2: "; PrintSized(i, r);  
+
+    if( !n_IsZero( i, r) )
+    {
+  #ifdef HAVE_RINGS
+      TS_ASSERT( n_DivBy(s, i, r) );
+  #endif
+       
+      res = n_Div(s, i, r);
+    
+      clog<< "N*(N+1)/2: "; PrintSized(res, r);
+
+
+      number d = n_Sub(res, sum1, r);
+      TS_ASSERT( n_IsZeroDivisor(d, r) );
+      n_Delete(&d, r);
+      
+      if( n_GetChar(r) == 0 )
+      {
+        TS_ASSERT( n_Equal(sum1, res, r) );
+        TS_ASSERT( n_Equal(res, sum1, r) );
+      }
+    } else
+      TS_ASSERT_EQUALS( n_GetChar(r), 2);
+    
+
+    n_Delete(&s, r);  n_Delete(&i, r);
+    
+    n_Delete(&sum1, r); n_Delete(&res, r);    
+    
+
+    s = n_Init(0  , r);
+    ss = n_Init(0 , r);
+    for( int k = N; k >= 0; k-- )
+    {
+      i = n_Init(k, r);
+      ndInpAdd(s, i, r); // s += i
+
+      i = n_Neg(i, r);
+      ndInpAdd(ss, i, r); // ss -= i
+      
+      n_Delete(&i, r);    
+    }
+    clog<< "ss: "; PrintSized(ss, r);  
+
+    ss = n_Neg(ss, r); // ss = -ss
+   
+    clog<< "real sum    : "; PrintSized(s, r);
+    clog<< "real sum(--): "; PrintSized(ss, r);  
+
+    TS_ASSERT( n_Equal(s, ss, r) );
+    TS_ASSERT( n_Equal(ss, s, r) );
+
+    n_Delete(&s, r);    
+    n_Delete(&ss, r);    
+
+    clog << ( " >>> TEST DONE!" );
+    clog << endl;
+
+  }
 public:
   void test_Z13_t()
@@ -388 +557 @@
 
     TS_ASSERT_DIFFERS( cf->cfCoeffWrite, NULLp );
-    
-    
   
     if( cf->cfCoeffWrite != NULL )
@@ -396 +563 @@
       n_CoeffWrite(cf); PrintLn();
     }
+    
+    // some tests for the coefficient field represented by cf:
+    TestArithCf(cf);
+    TestSumCf(cf, 10);
+    TestSumCf(cf, 100);
+    TestSumCf(cf, 101);
+    TestSumCf(cf, 1001);
+    TestSumCf(cf, 9000);
 
     clog << "Finally create the polynomial ring (Q[a]/<a2+1>)[x, y]..."
@@ -416 +591 @@
     TS_ASSERT_EQUALS(rVar(s), 2);
 
-    //Test(s);
+    Test(s);
 
     rDelete(s); // kills 'cf' and 'r' as well