Changeset 513673 in git for Singular/LIB/realclassify.lib

Timestamp:

Oct 15, 2013, 3:15:25 PM (11 years ago)

Author:

Andreas Steenpass <steenpass@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

3630419158783d9329a9e0d3d501208a60c71b10

Parents:

c21a6b6ea8d600f0dee5b13cbb0c042d2fe7e574

git-author:

Andreas Steenpass <steenpass@mathematik.uni-kl.de>2013-10-15 15:15:25+02:00

git-committer:

Andreas Steenpass <steenpass@mathematik.uni-kl.de>2013-10-15 15:19:17+02:00

Message:

chg: update realclassify.lib

The implementation of the splitting lemma has been changed.

File:

• Singular/LIB/realclassify.lib (modified) (5 diffs)

Singular/LIB/realclassify.lib

@@ -28 +28 @@
 ";
 
+LIB "linalg.lib";
 LIB "elim.lib";
 LIB "primdec.lib";
@@ -475 +476 @@
 USAGE:    realmorsesplit(f[, mu]); f poly, mu int
 RETURN:   a list consisting of the corank of f, the inertia index, an upper
-          bound for the determinacy, the residual form of f and
-          the transformation
+          bound for the determinacy, and the residual form of f
 NOTE:     The characteristic of the basering must be zero, the monomial order
           must be local, f must be contained in maxideal(2) and the Milnor
@@ -486 +486 @@
 EXAMPLE:  example morsesplit; shows an example"
 {
-  /* auxiliary variables */
-  int i, j;
+  int i;
 
   /* error check */
@@ -532 +531 @@
   }
 
-  /* preliminary stuff */
-  list S;
+  /* compute the determinacy */
   int k = determinacy(f, mu);
   f = jet(f, k);
-  def br = basering;
-  map Phi = br, maxideal(1);
-  map phi;
-  poly a, p, r;
-
-  /* treat the variables one by one */
+
+  /* get jet(f, 2) right */
+  matrix H = concat(jet(jacob(jacob(f)), 0)/2, unitmat(n));
+  H = sym_reduce(H);
+  intvec perm_zero;
+  intvec perm_neg;
+  intvec perm_pos;
+  int c;
+  int lambda;
   for(i = 1; i <= n; i++)
   {
-    if(jet(f, 2)/var(i) == 0)
-    {
-      S = insert(S, i);
-    }
-    else
-    {
-      f, a, p, r = rewriteformorsesplit(f, k, i);
-      if(jet(a, 0) == 0)
-      {
-        for(j = i+1; j <= n; j++)
-        {
-          if(jet(f, 2)/(var(i)*var(j)) != 0)
-          {
-            break;
-          }
-        }
-        phi = br, maxideal(1);
-        phi[j] = var(j)+var(i);
-        Phi = phi(Phi);
-        f = phi(f);
-      }
-      f, a, p, r = rewriteformorsesplit(f, k, i);
-      while(p != 0)
-      {
-        phi = br, maxideal(1);
-        phi[i] = var(i)-p/(2*jet(a, 0));
-        Phi = phi(Phi);
-        f = phi(f);
-        f, a, p, r = rewriteformorsesplit(f, k, i);
-      }
-    }
-  }
-
-  /* sort variables according to corank */
-  int cr = size(S);
-  phi = br, 0:n;
-  j = 1;
-  for(i = size(S); i > 0; i--)
-  {
-    phi[S[i]] = var(j);
-    j++;
-  }
+    if(H[i, i] == 0)
+    {
+      perm_zero = perm_zero, i;
+      c++;
+    }
+    if(H[i, i] < 0)
+    {
+      perm_neg = perm_neg, i;
+      lambda++;
+    }
+    if(H[i, i] > 0)
+    {
+      perm_pos = perm_pos, i;
+    }
+  }
+  intvec perm;
+  if(size(perm_zero) > 1)
+  {
+    perm = perm, perm_zero[2..size(perm_zero)];
+  }
+  if(size(perm_neg) > 1)
+  {
+    perm = perm, perm_neg[2..size(perm_neg)];
+  }
+  if(size(perm_pos) > 1)
+  {
+    perm = perm, perm_pos[2..size(perm_pos)];
+  }
+  perm = perm[2..size(perm)];
+  matrix T[n][n];
+  matrix D[1][n];
   for(i = 1; i <= n; i++)
   {
-    if(phi[i] == 0)
-    {
-      phi[i] = var(j);
-      j++;
-    }
-  }
-  Phi = phi(Phi);
+    T[1..n, i] = H[perm[i], (n+1)..(2*n)];
+    D[1, i] = H[perm[i], perm[i]];
+  }
+  map phi = basering, matrix(maxideal(1))*transpose(T);
   f = phi(f);
-
-  /* compute the inertia index lambda */
-  int lambda;
-  list negCoeff, posCoeff;
-  number ai;
-  poly f2 = jet(f, 2);
-  for(i = 1; i <= n; i++)
-  {
-    ai = number(f2/var(i)^2);
-    if(ai < 0)
-    {
-      lambda++;
-      negCoeff = insert(negCoeff, i);
-    }
-    if(ai > 0)
-    {
-      posCoeff = insert(posCoeff, i);
-    }
-  }
-
-  /* sort variables according to lambda */
-  phi = br, maxideal(1);
-  j = cr+1;
-  for(i = size(negCoeff); i > 0; i--)
-  {
-    phi[negCoeff[i]] = var(j);
-    j++;
-  }
-  for(i = size(posCoeff); i > 0; i--)
-  {
-    phi[posCoeff[i]] = var(j);
-    j++;
-  }
-  Phi = phi(Phi);
-  f = phi(f);
-
-  /* compute residual form */
-  phi = br, maxideal(1);
-  for(i = size(S)+1; i <= n; i++)
-  {
-    phi[i] = 0;
-  }
-  f = phi(f);
-
-  return(list(cr, lambda, k, f, Phi));
+  f = jet(f, k);
+
+  /* separate the variables */
+  phi = basering, maxideal(1);
+  map corank_part = basering, maxideal(1);
+  for(i = c+1; i <= n; i++)
+  {
+    corank_part[i] = 0;
+  }
+  poly h = f-jet(f, 2)-corank_part(f);
+  poly hi;
+  while(h != 0)
+  {
+    for(i = c+1; i <= n; i++)
+    {
+      hi = h/var(i);
+      phi[i] = var(i)-hi/(2*D[1, i]);
+      h = h-hi*var(i);
+    }
+    f = phi(f);
+    f = jet(f, k);
+    h = f-jet(f, 2)-corank_part(f);
+  }
+  poly g = cleardenom(f-jet(f, 2));
+
+  return(list(c, lambda, k, g));
 }
 example
@@ -649 +617 @@
   poly f = (x2+3y-2z)^2+xyz-(x-y3+x2z3)^3;
   realmorsesplit(f);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+/*
+   symmetric Gauss algorithm
+
+   If A is not a square matrix, then the largest upper or left submatrix
+   is assumed to be symmetric.
+*/
+proc sym_reduce(matrix A)
+{
+  int r = nrows(A);
+  int c = ncols(A);
+  int n = r;
+  if(n > c)
+  {
+    n = c;
+  }
+  poly q;
+  int i, j;
+  for(i = 1; i <= n; i++)
+  {
+    for(j = i+1; j <= n; j++)
+    {
+      if(A[i, j] != 0)
+      {
+        while(A[i, i] == 0)
+        {
+          A[1..r, i] = A[1..r, i]+A[1..r, j];
+          A[i, 1..c] = A[i, 1..c]+A[j, 1..c];
+        }
+        q = A[i, j]/A[i, i];
+        A[1..r, j] = A[1..r, j]-q*A[1..r, i];
+        A[j, 1..c] = A[j, 1..c]-q*A[i, 1..c];
+      }
+    }
+  }
+  return(A);
 }