Changeset 47b559d in git


Ignore:
Timestamp:
Dec 2, 2010, 12:37:29 PM (13 years ago)
Author:
Frank Seelisch <seelisch@…>
Branches:
(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
Children:
2f2562016613663d129949817c3520b92f72c095
Parents:
35715cbb715ef8f8977fca4b923054147968c0d1
Message:
finished hensel factorization (provided as sys command in extra.cc)

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

Legend:

Unmodified
Added
Removed

  • Singular/extra.cc

    r35715c r47b559d  
    633633      if(strcmp(sys_cmd, "henselfactors")==0)
    634634      {
    635         if ((h != NULL) && (h->Typ() == POLY_CMD) &&
    636             (h->next != NULL) && (h->next->Typ() == POLY_CMD) &&
     635        if ((h != NULL) && (h->Typ() == INT_CMD) &&
     636            (h->next != NULL) && (h->next->Typ() == INT_CMD) &&
    637637            (h->next->next != NULL) && (h->next->next->Typ() == POLY_CMD) &&
    638638            (h->next->next->next != NULL) &&
    639                                      (h->next->next->next->Typ() == INT_CMD))
    640         {
    641           poly hh = (poly)h->Data();
    642           poly f0 = (poly)h->next->Data();
    643           poly g0 = (poly)h->next->next->Data();
    644           int d   = (int)(long)h->next->next->next->Data();
     639               (h->next->next->next->Typ() == POLY_CMD) &&
     640            (h->next->next->next->next != NULL) &&
     641               (h->next->next->next->next->Typ() == POLY_CMD) &&
     642            (h->next->next->next->next->next != NULL) &&
     643               (h->next->next->next->next->next->Typ() == INT_CMD) &&
     644            (h->next->next->next->next->next->next == NULL))
     645        {
     646          int xIndex = (int)(long)h->Data();
     647          int yIndex = (int)(long)h->next->Data();
     648          poly hh    = (poly)h->next->next->Data();
     649          poly f0    = (poly)h->next->next->next->Data();
     650          poly g0    = (poly)h->next->next->next->next->Data();
     651          int d      = (int)(long)h->next->next->next->next->next->Data();
    645652          poly f; poly g;
    646           henselFactors(hh, f0, g0, d, f, g);
     653          henselFactors(xIndex, yIndex, hh, f0, g0, d, f, g);
    647654          lists L = (lists)omAllocBin(slists_bin);
    648655          L->Init(2);
    … …  
    655662        else
    656663        {
    657           Werror( "expected argument list (poly, poly, poly, int)");
     664          Werror( "expected argument list (int, int, poly, poly, poly, int)");
    658665          return TRUE;
    659666        }

  • kernel/linearAlgebra.cc

    r35715c r47b559d  
    12631263}
    12641264
    1265 /* This codes makes no assumption about the monomial ordering in the current
    1266    base ring. */
    1267 void henselFactors(const poly h, const poly f0, const poly g0, const int d,
    1268                    poly &f, poly &g)
     1265/* This codes assumes that there are at least two variables in the current
     1266   base ring. No assumption is made regarding the monomial ordering. */
     1267void henselFactors(const int xIndex, const int yIndex, const poly h,
     1268                   const poly f0, const poly g0, const int d, poly &f, poly &g)
    12691269{
    12701270  int n = (int)pDeg(f0);
    12711271  int m = (int)pDeg(g0);
    1272   matrix fMat = mpNew(n, d + 1); /* fMat[i, j] is the coefficient of
    1273                                     x^(j-1) * y^(i-1) in f */
    1274   matrix gMat = mpNew(m, d + 1); /* gMat[i, j] is the coefficient of
    1275                                     x^(j-1) * y^(i-1) in g */
     1272  matrix fMat = mpNew(n + 1, d + 1);   /* fMat[i, j] is the coefficient of
     1273                                          x^(j-1) * y^(i-1) in f */
     1274  matrix gMat = mpNew(m + 1, d + 1);   /* gMat[i, j] is the coefficient of
     1275                                          x^(j-1) * y^(i-1) in g */
     1276  matrix hMat = mpNew(n + m, d + 1);   /* hMat[i, j] is the coefficient of
     1277                                          x^(j-1) * y^(i-1) in h */
    12761278 
    1277   f = pOne(); pSetExp(f, 2, n); pSetm(f);   // initialization: f = y^n
    1278   g = pOne(); pSetExp(g, 2, m); pSetm(g);   // initialization: g = y^m
    1279  
    1280   /* initial step: read off coefficients of f0 and g0 */
    1281   poly l = pAdd(pCopy(f0), pNeg(pCopy(f)));   // l = f0 - f = f0 - y^n
     1279  /* initial step: read off coefficients of f0, g0, and h */
     1280  poly l = f0;
    12821281  while (l != NULL)
    12831282  {
    1284     int e = pGetExp(l, 2);
     1283    int e = pGetExp(l, yIndex);
    12851284    number c = pGetCoeff(l);
    12861285    poly matEntry = pOne(); pSetCoeff(matEntry, nCopy(c));
    … …  
    12881287    l = pNext(l);
    12891288  }
    1290   l = pAdd(pCopy(g0), pNeg(pCopy(g)));        // l = g0 - g = g0 - y^m
     1289  l = g0;
    12911290  while (l != NULL)
    12921291  {
    1293     int e = pGetExp(l, 2);
     1292    int e = pGetExp(l, yIndex);
    12941293    number c = pGetCoeff(l);
    12951294    poly matEntry = pOne(); pSetCoeff(matEntry, nCopy(c));
    … …  
    12971296    l = pNext(l);
    12981297  }
     1298  l = h;
     1299  while (l != NULL)
     1300  {
     1301    int eX = pGetExp(l, xIndex);
     1302    if (eX <= d) /* we forget about terms with x-exponents > d */
     1303    {
     1304      int eY = pGetExp(l, yIndex);
     1305      if (eY <= n + m - 1) /* we forget about the leading term y^(n + m) */
     1306      {
     1307        number c = pGetCoeff(l);
     1308        poly matEntry = pOne(); pSetCoeff(matEntry, nCopy(c));
     1309        MATELEM(hMat, eY + 1, eX + 1) = matEntry;
     1310      }
     1311    }
     1312    l = pNext(l);
     1313  }
    12991314 
    13001315  /* loop for computing entries of fMat and gMat in columns 2..(d+1)
    13011316     which correspond to the x-exponents 1..d */
    1302   int p = m; if (n < m) p = n;   // p = min{m, n}
    1303   for (int exp = 1; exp <= d; exp++)
    1304   {
    1305     matrix aMat = mpNew(p, n + m);   // container for matrix A of linear system
    1306     matrix xVec = mpNew(n + m, 1);   // container for unknowns of linear system
    1307     matrix bVec = mpNew(p, 1);       // container for result vector
     1317  for (int xExp = 1; xExp <= d; xExp++)
     1318  {
     1319    /* We are going to setup and solve a linear equation system A * x = b,
     1320       where the resulting x will give us all matrix entries fMat[i, xExp],
     1321       i = 1, 2, ..., n; and gMat[j, xExp], j = 1, 2, ..., m.
     1322       (Theory ensures that the system is uniquely solvable.) */
     1323    matrix aMat = mpNew(n + m, n + m); /* A */
     1324    matrix xVec = mpNew(n + m, 1);     /* x */
     1325    matrix bVec = mpNew(n + m, 1);     /* b */
    13081326   
    1309     // continue here
    1310    
     1327    /* setup A */
     1328    for (int row = 1; row <= n + m; row++)
     1329    {
     1330      int k = row;
     1331      for (int col = 1; col <= n; col++)
     1332      {
     1333        if (k <= m + 1) MATELEM(aMat, row, col) = pCopy(MATELEM(gMat, k, 1));
     1334        k--;
     1335        if (k == 0) break;
     1336      }
     1337      k = row;
     1338      for (int col = n + 1; col <= n + m; col++)
     1339      {
     1340        if (k <= n + 1) MATELEM(aMat, row, col) = pCopy(MATELEM(fMat, k, 1));
     1341        k--;
     1342        if (k == 0) break;
     1343      }
     1344    }
     1345
     1346    /* setup b */
     1347    for (int row = 1; row <= n + m; row++)
     1348    {
     1349      poly p = pCopy(MATELEM(hMat, row, xExp + 1));
     1350      for (int j = 1; j <= row; j++)
     1351      {
     1352        if (j > m + 1) break;
     1353        if (row + 1 - j <= n + 1)
     1354        {
     1355          for (int i = 1; i <= xExp + 1; i++)
     1356          {
     1357            poly q = ppMult_qq(MATELEM(fMat, row + 1 - j, xExp + 2 - i),
     1358                               MATELEM(gMat, j, i));
     1359            q = pNeg(q);
     1360            p = pAdd(p, q);
     1361          }
     1362        }
     1363      }
     1364      MATELEM(bVec, row, 1) = p;
     1365    }
     1366
     1367    /* computation of x */
     1368    matrix pMat; matrix lMat; matrix uMat; matrix wMat;
     1369    luDecomp(aMat, pMat, lMat, uMat);
     1370    luSolveViaLUDecomp(pMat, lMat, uMat, bVec, xVec, wMat);
     1371
     1372    /* fill (xExp + 1)-columns of fMat and gMat */
     1373    for (int row = 1; row <= n; row++)
     1374      MATELEM(fMat, row, xExp + 1) = pCopy(MATELEM(xVec, row, 1));
     1375    for (int row = 1; row <= m; row++)
     1376      MATELEM(gMat, row, xExp + 1) = pCopy(MATELEM(xVec, n + row, 1));
    13111377   
    13121378    /* clean-up */
    1313     idDelete((ideal*)&aMat);
    1314     idDelete((ideal*)&xVec);
    1315     idDelete((ideal*)&bVec);
     1379    idDelete((ideal*)&aMat); idDelete((ideal*)&xVec); idDelete((ideal*)&bVec);
     1380    idDelete((ideal*)&pMat); idDelete((ideal*)&lMat); idDelete((ideal*)&uMat);
     1381    idDelete((ideal*)&wMat);
    13161382  }
    13171383 
    13181384  /* build f and g from fMat and gMat, respectively */
    1319   for (int i = 0; i < n; i++)
     1385  f = NULL;
     1386  for (int i = 0; i <= n; i++)
    13201387    for (int j = 0; j <= d; j++)
    13211388    {
    … …  
    13231390      if (term != NULL)
    13241391      {
    1325         pSetExp(term, 1, j);                            // ... * x^j
    1326         pSetExp(term, 2, i);                            // ... * y^i
     1392        pSetExp(term, xIndex, j);                       // ... * x^j
     1393        pSetExp(term, yIndex, i);                       // ... * y^i
    13271394        pSetm(term);
    13281395        f = pAdd(f, term);
    13291396      }
    13301397    }
    1331   for (int i = 0; i < m; i++)
     1398  g = NULL;
     1399  for (int i = 0; i <= m; i++)
    13321400    for (int j = 0; j <= d; j++)
    13331401    {
    … …  
    13351403      if (term != NULL)
    13361404      {
    1337         pSetExp(term, 1, j);                            // ... * x^j
    1338         pSetExp(term, 2, i);                            // ... * y^i
     1405        pSetExp(term, xIndex, j);                       // ... * x^j
     1406        pSetExp(term, yIndex, i);                       // ... * y^i
    13391407        pSetm(term);
    13401408        g = pAdd(g, term);
    … …  
    13451413  idDelete((ideal*)&fMat);
    13461414  idDelete((ideal*)&gMat);
    1347 }
    1348 
     1415  idDelete((ideal*)&hMat);
     1416}
     1417

  • kernel/linearAlgebra.h

    r35715c r47b559d  
    470470 * all terms of h with x-degree larger than d can be ignored due to (*).
    471471 * The method expects the ground ring to contain at least two variables; then
    472  * x is the first ring variable and y the second one.
     472 * x is the first ring variable, specified by the input index xIndex, and y the
     473 * second one, specified by yIndex.
    473474 *
    474475 * This code was place here since the algorithm works by successively solving
    … …  
    478479 **/
    479480void henselFactors(
     481       const int xIndex,  /**< [in]  index of x in {1, ..., nvars(basering)} */
     482       const int yIndex,  /**< [in]  index of y in {1, ..., nvars(basering)} */
    480483       const poly h,      /**< [in]  the polynomial h(x, y) to be factorized */
    481484       const poly f0,     /**< [in]  the first univariate factor of h(0, y)  */
Note: See TracChangeset for help on using the changeset viewer.