Changeset 35715c in git

Timestamp:

Dec 1, 2010, 5:01:09 PM (13 years ago)

Author:

Frank Seelisch <seelisch@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

47b559d8a480780a6380da801528adb53dca1000

Parents:

af71455603162123e6412ac27852bb08ea8f5126

Message:

further development of hensel factorization

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

Location:

kernel

Files:

• linearAlgebra.cc (modified) (2 diffs)
• linearAlgebra.h (modified) (3 diffs)

kernel/linearAlgebra.cc

@@ -861 +861 @@
   }
   nDelete(&denominator);
-                          /*
-                          // inspect the following for inaccuracies in the computation;
-                          // norm of u should be 2 by construction
-                                poly o = NULL;
-                                for (int j = 1; j <= rr; j++)
-                                  o = pAdd(o, ppMult_qq(MATELEM(uVec, j, 1), MATELEM(uVec, j, 1)));     
-                                printf("norm of u = %s\n", pString(o));
-                                printMatrix(vVec); printMatrix(uVec);
-                                pDelete(&o);
-                          */
+
   /* finished building vector u; now turn to pMat */
   pMat = mpNew(rr, rr);
@@ -1272 +1263 @@
 }
 
+/* This codes makes no assumption about the monomial ordering in the current
+   base ring. */
 void henselFactors(const poly h, const poly f0, const poly g0, const int d,
                    poly &f, poly &g)
 {
-  f = NULL;
-  g = NULL;
-}
-
+  int n = (int)pDeg(f0);
+  int m = (int)pDeg(g0);
+  matrix fMat = mpNew(n, d + 1); /* fMat[i, j] is the coefficient of
+                                    x^(j-1) * y^(i-1) in f */
+  matrix gMat = mpNew(m, d + 1); /* gMat[i, j] is the coefficient of
+                                    x^(j-1) * y^(i-1) in g */
+  
+  f = pOne(); pSetExp(f, 2, n); pSetm(f);   // initialization: f = y^n
+  g = pOne(); pSetExp(g, 2, m); pSetm(g);   // initialization: g = y^m
+  
+  /* initial step: read off coefficients of f0 and g0 */
+  poly l = pAdd(pCopy(f0), pNeg(pCopy(f)));   // l = f0 - f = f0 - y^n
+  while (l != NULL)
+  {
+    int e = pGetExp(l, 2);
+    number c = pGetCoeff(l);
+    poly matEntry = pOne(); pSetCoeff(matEntry, nCopy(c));
+    MATELEM(fMat, e + 1, 1) = matEntry;
+    l = pNext(l);
+  }
+  l = pAdd(pCopy(g0), pNeg(pCopy(g)));        // l = g0 - g = g0 - y^m
+  while (l != NULL)
+  {
+    int e = pGetExp(l, 2);
+    number c = pGetCoeff(l);
+    poly matEntry = pOne(); pSetCoeff(matEntry, nCopy(c));
+    MATELEM(gMat, e + 1, 1) = matEntry;
+    l = pNext(l);
+  }
+  
+  /* loop for computing entries of fMat and gMat in columns 2..(d+1)
+     which correspond to the x-exponents 1..d */
+  int p = m; if (n < m) p = n;   // p = min{m, n}
+  for (int exp = 1; exp <= d; exp++)
+  {
+    matrix aMat = mpNew(p, n + m);   // container for matrix A of linear system
+    matrix xVec = mpNew(n + m, 1);   // container for unknowns of linear system
+    matrix bVec = mpNew(p, 1);       // container for result vector
+    
+    // continue here
+    
+    
+    /* clean-up */
+    idDelete((ideal*)&aMat);
+    idDelete((ideal*)&xVec);
+    idDelete((ideal*)&bVec);
+  }
+ 
+  /* build f and g from fMat and gMat, respectively */
+  for (int i = 0; i < n; i++)
+    for (int j = 0; j <= d; j++)
+    {
+      poly term = pCopy(MATELEM(fMat, i + 1, j + 1));   // coeff
+      if (term != NULL)
+      {
+        pSetExp(term, 1, j);                            // ... * x^j
+        pSetExp(term, 2, i);                            // ... * y^i
+        pSetm(term);
+        f = pAdd(f, term);
+      }
+    }
+  for (int i = 0; i < m; i++)
+    for (int j = 0; j <= d; j++)
+    {
+      poly term = pCopy(MATELEM(gMat, i + 1, j + 1));   // coeff
+      if (term != NULL)
+      {
+        pSetExp(term, 1, j);                            // ... * x^j
+        pSetExp(term, 2, i);                            // ... * y^i
+        pSetm(term);
+        g = pAdd(g, term);
+      }
+    }
+  
+  /* clean-up */
+  idDelete((ideal*)&fMat);
+  idDelete((ideal*)&gMat);
+}
+

kernel/linearAlgebra.h

@@ -219 +219 @@
 
 /**
- * Solves the linear system A*x = b, where A is an (n x n)-matrix
+ * Solves the linear system A * x = b, where A is an (m x n)-matrix
  * which is given by its LU-decomposition.
  *
@@ -461 +461 @@
  * polynomial in y of degree m + n with coefficients in K[[x]]. Suppose there
  * are two monic factors f_0(y) (of degree n) and g_0(y) of degree (m) such
- * that h(0, y) = f_0(y) * g_0(y), and an integer d >= 0. Then there are
- * monic polynomials in y with coefficients in K[[x]], namely f(x, y) of
- * degree n and g(x, y) of degree m such that
+ * that h(0, y) = f_0(y) * g_0(y) and <f_0, g_0> = K[y]. Fix an integer d >= 0.
+ * Then there are monic polynomials in y with coefficients in K[[x]], namely
+ * f(x, y) of degree n and g(x, y) of degree m such that
  *    h(x, y) = f(x, y) * g(x, y) modulo <x^(d+1)>   (*).
  *
@@ -471 +471 @@
  * The method expects the ground ring to contain at least two variables; then
  * x is the first ring variable and y the second one.
+ *
+ * This code was place here since the algorithm works by successively solving
+ * d linear equation systems. It is hence an application of other methods
+ * defined in this h-file and its corresponding cc-file.
+ * 
  **/
 void henselFactors(