Changeset 139b67 in git


Ignore:
Timestamp:
Sep 13, 2010, 4:37:37 PM (13 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
Children:
231b0fd0173552d615fba69018a64ff5510211be
Parents:
ce35eb186c60fc2922d2d106228a4d3a2cb87a5b
Message:
new routines for two factor hensel lifting and some bug fixes


git-svn-id: file:///usr/local/Singular/svn/trunk@13183 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
factory
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • factory/cf_gcd_smallp.cc

    rce35eb r139b67  
    904904  bool inextensionextension= false;
    905905  top_level= false;
    906   CanonicalForm CF_buf;
    907906  CFList source, dest;
    908   CanonicalForm gcdcheck;
    909907  do
    910908  {
    … …  
    942940      CFList list;
    943941      CanonicalForm mipo;
    944       int deg= 3;
     942      int deg= 2;
    945943      do {
    946944        mipo= randomIrredpoly (deg, x);

  • factory/cf_map_ext.cc

    rce35eb r139b67  
    340340      return 0;
    341341  } while (1);
    342   zz_pE::init (NTL_mipo);
    343   zz_pEX NTL_alpha_mipo= convertFacCF2NTLzz_pEX (mipo, NTL_mipo);
    344   zz_pE root= FindRoot (NTL_alpha_mipo);
     342 
     343  zz_pX alpha_mipo= convertFacCF2NTLzzpX (mipo);
     344  zz_pE::init (alpha_mipo);
     345  zz_pEX NTL_beta_mipo= to_zz_pEX (NTL_mipo);
     346  zz_pE root= FindRoot (NTL_beta_mipo);
    345347  return convertNTLzzpE2CF (root, alpha);
    346348}

  • factory/facHensel.cc

    rce35eb r139b67  
    2626#include "facHensel.h"
    2727#include "cf_util.h"
     28#include "fac_util.h"
     29#include "cf_algorithm.h"
    2830
    2931#ifdef HAVE_NTL
    3032#include <NTL/lzz_pEX.h>
    3133#include "NTLconvert.h"
     34
     35CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c);
    3236
    3337static inline
    … …  
    920924void
    921925henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
    922               CFList& diophant, CFMatrix& M)
    923 {
    924   sortList (factors, Variable (1));
     926              CFList& diophant, CFMatrix& M, bool sort)
     927{
     928  if (sort)
     929    sortList (factors, Variable (1));
    925930  Pi= CFArray (factors.length() - 1);
    926931  CFListIterator j= factors;
    … …  
    933938  if (j.hasItem())
    934939    j++;
    935   for (j; j.hasItem(); j++, i++)
     940  for (; j.hasItem(); j++, i++)
    936941  {
    937942    Pi [i]= mulNTL (Pi [i - 1], j.getItem());
    … …  
    954959    k.getItem()= bufFactors[i];
    955960  factors.removeFirst();
     961  return;
     962}
     963
     964void
     965henselStep122 (const CanonicalForm& F, const CFList& factors,
     966              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
     967              CFArray& Pi, int j, const CFArray& LCs)
     968{
     969  Variable x= F.mvar();
     970  CanonicalForm xToJ= power (x, j);
     971 
     972  CanonicalForm E;
     973  // compute the error
     974  if (degree (Pi [factors.length() - 2], x) > 0)
     975    E= F[j] - Pi [factors.length() - 2] [j];
     976  else
     977    E= F[j];
     978
     979  CFArray buf= CFArray (diophant.length());
     980
     981  int k= 0;
     982  CanonicalForm remainder;
     983  // actual lifting
     984  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
     985  {
     986    if (degree (bufFactors[k], x) > 0)
     987      remainder= modNTL (E, bufFactors[k] [0]);
     988    else
     989      remainder= modNTL (E, bufFactors[k]);
     990    buf[k]= mulNTL (i.getItem(), remainder);
     991    if (degree (bufFactors[k], x) > 0)
     992      buf[k]= modNTL (buf[k], bufFactors[k] [0]);
     993    else
     994      buf[k]= modNTL (buf[k], bufFactors[k]);
     995  }
     996 
     997  for (k= 0; k < factors.length(); k++)
     998    bufFactors[k] += xToJ*buf[k];
     999   
     1000  // update Pi [0]
     1001  int degBuf0= degree (bufFactors[0], x);
     1002  int degBuf1= degree (bufFactors[1], x);
     1003  if (degBuf0 > 0 && degBuf1 > 0)
     1004  {
     1005    M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
     1006    if (j + 2 <= M.rows())
     1007      M (j + 2, 1)= mulNTL (bufFactors[0] [j + 1], bufFactors[1] [j + 1]);
     1008  }
     1009  CanonicalForm uIZeroJ;
     1010  if (degBuf0 > 0 && degBuf1 > 0)
     1011    uIZeroJ= mulNTL (bufFactors[0] [0], buf [1]) + mulNTL (bufFactors[1] [0], buf[0]);
     1012  else if (degBuf0 > 0)
     1013    uIZeroJ= mulNTL (buf[0], bufFactors[1]);
     1014  else if (degBuf1 > 0)
     1015    uIZeroJ= mulNTL (bufFactors[0], buf [1]);
     1016  else
     1017    uIZeroJ= 0;
     1018  Pi [0] += xToJ*uIZeroJ;
     1019
     1020  CFArray tmp= CFArray (factors.length() - 1);
     1021  for (k= 0; k < factors.length() - 1; k++)
     1022    tmp[k]= 0;
     1023  CFIterator one, two;
     1024  one= bufFactors [0];
     1025  two= bufFactors [1];
     1026  if (degBuf0 > 0 && degBuf1 > 0)
     1027  {
     1028    while (one.exp() > j) one++;
     1029    while (two.exp() > j) two++;
     1030    for (k= 1; k <= (int) ceil (j/2.0); k++)
     1031    { 
     1032      if (k != j - k + 1)
     1033      {
     1034        if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
     1035        {
     1036          tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), (bufFactors[1] [k] +
     1037                     two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
     1038          one++;
     1039          two++;
     1040        }
     1041        else if (one.exp() == j - k + 1)
     1042        {
     1043          tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), bufFactors[1] [k]) -
     1044                     M (k + 1, 1);
     1045          one++;
     1046        }
     1047        else if (two.exp() == j - k + 1)
     1048        {
     1049          tmp[0] += mulNTL (bufFactors[0] [k], (bufFactors[1] [k] + two.coeff())) -
     1050                    M (k + 1, 1);
     1051          two++;
     1052        }
     1053      }
     1054      else
     1055        tmp[0] += M (k + 1, 1);
     1056    }
     1057
     1058    if (j + 2 <= M.rows())
     1059    {
     1060      if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
     1061        tmp [0] += mulNTL ((bufFactors [0] [j + 1]+ bufFactors [0] [0]), (bufFactors [1] [j + 1] + bufFactors [1] [0]))
     1062               - M(1,1) - M (j + 2,1);
     1063      else if (degBuf0 >= j + 1)
     1064        tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1] [0]);
     1065      else if (degBuf1 >= j + 1)
     1066        tmp[0] += mulNTL (bufFactors [0] [0], bufFactors [1] [j + 1]);
     1067    }
     1068  }
     1069  Pi [0] += tmp[0]*xToJ*F.mvar();
     1070 
     1071  /*// update Pi [l]
     1072  int degPi, degBuf;
     1073  for (int l= 1; l < factors.length() - 1; l++)
     1074  {
     1075    degPi= degree (Pi [l - 1], x);
     1076    degBuf= degree (bufFactors[l + 1], x);
     1077    if (degPi > 0 && degBuf > 0)
     1078      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]); 
     1079    if (j == 1)
     1080    {
     1081      if (degPi > 0 && degBuf > 0)
     1082        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
     1083                  bufFactors[l + 1] [0] + buf[l + 1]) - M (j + 1, l +1) -
     1084                  M (1, l + 1));
     1085      else if (degPi > 0)
     1086        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1]));
     1087      else if (degBuf > 0)
     1088        Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1]));
     1089    }
     1090    else
     1091    {
     1092      if (degPi > 0 && degBuf > 0)
     1093      {
     1094        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
     1095        uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1]);
     1096      }
     1097      else if (degPi > 0)
     1098        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1]);
     1099      else if (degBuf > 0)
     1100      {
     1101        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
     1102        uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1]);
     1103      }
     1104      Pi[l] += xToJ*uIZeroJ;
     1105    }
     1106    one= bufFactors [l + 1];
     1107    two= Pi [l - 1];
     1108    if (two.exp() == j + 1)
     1109    {
     1110      if (degBuf > 0 && degPi > 0)
     1111      {
     1112          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0]);
     1113          two++;
     1114      }
     1115      else if (degPi > 0)
     1116      {     
     1117          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1]);
     1118          two++;
     1119      }
     1120    }
     1121    if (degBuf > 0 && degPi > 0)
     1122    {
     1123      for (k= 1; k <= (int) ceil (j/2.0); k++)
     1124      { 
     1125        if (k != j - k + 1)
     1126        {
     1127          if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
     1128          {
     1129            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), (Pi[l - 1] [k] +
     1130                       two.coeff())) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
     1131            one++;
     1132            two++;
     1133          }
     1134          else if (one.exp() == j - k + 1)
     1135          {
     1136            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), Pi[l - 1] [k]) -
     1137                       M (k + 1, l + 1);     
     1138            one++;
     1139          }
     1140          else if (two.exp() == j - k + 1)
     1141          {
     1142            tmp[l] += mulNTL (bufFactors[l + 1] [k], (Pi[l - 1] [k] + two.coeff())) -
     1143                      M (k + 1, l + 1);
     1144            two++;
     1145          }
     1146        }
     1147        else
     1148          tmp[l] += M (k + 1, l + 1);
     1149      } 
     1150    } 
     1151    Pi[l] += tmp[l]*xToJ*F.mvar();
     1152  }*/
     1153  return;
     1154}
     1155
     1156void
     1157henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
     1158              CFList& diophant, CFMatrix& M, const CFArray& LCs, bool sort)
     1159{
     1160  if (sort)
     1161    sortList (factors, Variable (1));
     1162  Pi= CFArray (factors.length() - 2);
     1163  CFList bufFactors2= factors;
     1164  bufFactors2.removeFirst();
     1165  diophant.removeFirst();
     1166  CFListIterator iter= diophant;
     1167  CanonicalForm test= 0, test2= 0;
     1168  CanonicalForm s,t;
     1169  extgcd (bufFactors2.getFirst(), bufFactors2.getLast(), s, t);
     1170  diophant= CFList();
     1171  diophant.append (t);
     1172  diophant.append (s);
     1173  DEBOUTLN (cerr, "diophant= " << diophant);
     1174
     1175  CFArray bufFactors= CFArray (bufFactors2.length());
     1176  int i= 0;
     1177  for (CFListIterator k= bufFactors2; k.hasItem(); i++, k++)
     1178    bufFactors[i]= replaceLc (k.getItem(), LCs [i]);
     1179
     1180  Variable x= F.mvar();
     1181  if (degree (bufFactors[0], x) > 0 && degree (bufFactors [1], x) > 0)
     1182  {
     1183    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1] [0]);
     1184    Pi [0]= M (1, 1) + (mulNTL (bufFactors [0] [1], bufFactors[1] [0]) +
     1185                        mulNTL (bufFactors [0] [0], bufFactors [1] [1]))*x;
     1186  }
     1187  else if (degree (bufFactors[0], x) > 0)
     1188  {
     1189    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1]);
     1190    Pi [0]= M (1, 1) +
     1191            mulNTL (bufFactors [0] [1], bufFactors[1])*x;
     1192  }
     1193  else if (degree (bufFactors[1], x) > 0)
     1194  {
     1195    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1] [0]);
     1196    Pi [0]= M (1, 1) +
     1197            mulNTL (bufFactors [0], bufFactors[1] [1])*x;
     1198  }
     1199  else
     1200  {
     1201    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1]);
     1202    Pi [0]= M (1, 1);
     1203  }         
     1204
     1205  for (i= 1; i < l; i++)
     1206    henselStep122 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
     1207
     1208  factors= CFList();
     1209  for (i= 0; i < bufFactors.size(); i++)
     1210    factors.append (bufFactors[i]);
    9561211  return;
    9571212}
    … …  
    11331388  for (int i= 1; i < d; i++)
    11341389  {
     1390    CanonicalForm kk= 0;
    11351391    if (degree (e, y) > 0)
    11361392      coeffE= e[i];
    … …  
    11561412      }
    11571413    }
    1158  
     1414    if (e.isZero())
     1415      break;
     1416  }
     1417
     1418#ifdef DEBUGOUTPUT
     1419    CanonicalForm test= 0;
     1420    j= p;
     1421    for (CFListIterator i= result; i.hasItem(); i++, j++)
     1422      test += mod (i.getItem()*j.getItem(), power (y, d));
     1423    DEBOUTLN (cerr, "test= " << test);
     1424#endif
     1425  return result;
     1426}
     1427
     1428/// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$
     1429/// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$
     1430static inline
     1431CFList
     1432diophantine (const CFList& recResult, const CFList& factors,
     1433             const CFList& products, const CFList& M, const CanonicalForm& E)
     1434{
     1435  if (M.isEmpty())
     1436  {
     1437    CFList result;
     1438    CFListIterator j= factors;
     1439    CanonicalForm buf;
     1440    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
     1441    {
     1442      ASSERT (E.isUnivariate() || E.inCoeffDomain(),
     1443              "constant or univariate poly expected");
     1444      ASSERT (i.getItem().isUnivariate() || i.getItem().inCoeffDomain(),
     1445              "constant or univariate poly expected");
     1446      ASSERT (j.getItem().isUnivariate() || j.getItem().inCoeffDomain(),
     1447              "constant or univariate poly expected");
     1448      buf= mulNTL (E, i.getItem());
     1449      result.append (modNTL (buf, j.getItem()));
     1450    }
     1451    return result;
     1452  }
     1453  Variable y= M.getLast().mvar();
     1454  CFList bufFactors= factors;
     1455  for (CFListIterator i= bufFactors; i.hasItem(); i++)
     1456    i.getItem()= mod (i.getItem(), y);
     1457  CFList bufProducts= products;
     1458  for (CFListIterator i= bufProducts; i.hasItem(); i++)
     1459    i.getItem()= mod (i.getItem(), y);
     1460  CFList buf= M;
     1461  buf.removeLast();
     1462  CanonicalForm bufE= mod (E, y);
     1463  CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
     1464                                      bufE);
     1465
     1466  CanonicalForm e= E;
     1467  CFListIterator j= products;
     1468  for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++)
     1469    e -= mulMod (i.getItem(), j.getItem(), M);
     1470
     1471  CFList result= recDiophantine;
     1472  int d= degree (M.getLast());
     1473  CanonicalForm coeffE;
     1474  for (int i= 1; i < d; i++)
     1475  {
     1476    if (degree (e, y) > 0)
     1477      coeffE= e[i];
     1478    else
     1479      coeffE= 0;
     1480    if (!coeffE.isZero())
     1481    {
     1482      CFListIterator k= result;
     1483      recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
     1484                                   coeffE);
     1485      CFListIterator l= products;
     1486      for (j= recDiophantine; j.hasItem(); j++, k++, l++)
     1487      {
     1488        k.getItem() += j.getItem()*power (y, i);
     1489        e -= mulMod (j.getItem()*power (y, i), l.getItem(), M);
     1490      }
     1491    }
     1492    if (e.isZero())
     1493      break;
     1494  }
     1495  return result;
     1496}
     1497
     1498// non monic case
     1499static inline
     1500CFList
     1501multiRecDiophantine2 (const CanonicalForm& F, const CFList& factors,
     1502                     const CFList& recResult, const CFList& M, const int d)
     1503{
     1504  Variable y= F.mvar();
     1505  CFList result;
     1506  CFListIterator i;
     1507  CanonicalForm e= 1;
     1508  CFListIterator j= factors;
     1509  CFList p;
     1510  CFArray bufFactors= CFArray (factors.length());
     1511  CanonicalForm yToD= power (y, d);
     1512  int k= 0;
     1513  for (CFListIterator i= factors; i.hasItem(); i++, k++)
     1514    bufFactors [k]= i.getItem();
     1515  CanonicalForm b;
     1516  CFList buf= M;
     1517  buf.removeLast();
     1518  buf.append (yToD);
     1519  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
     1520  {
     1521    b= 1;
     1522    if (fdivides (bufFactors[k], F))
     1523      b= F/bufFactors[k];
     1524    else
     1525    {
     1526      for (int l= 0; l < factors.length(); l++)
     1527      {
     1528        if (l == k)
     1529          continue;
     1530        else
     1531        {
     1532          b= mulMod (b, bufFactors[l], buf);
     1533        }
     1534      }
     1535    }
     1536    p.append (b);
     1537  } 
     1538  j= p;
     1539  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
     1540    e -= mulMod (i.getItem(), j.getItem(), M);
     1541  if (e.isZero())
     1542    return recResult;
     1543  CanonicalForm coeffE;
     1544  result= recResult;
     1545  CanonicalForm g;
     1546  buf.removeLast();
     1547  for (int i= 1; i < d; i++)
     1548  {
     1549    if (degree (e, y) > 0)
     1550      coeffE= e[i];
     1551    else
     1552      coeffE= 0;
     1553    if (!coeffE.isZero())
     1554    {
     1555      CFListIterator k= result;
     1556      CFListIterator l= p;
     1557      j= recResult;
     1558      int ii= 0;
     1559      CanonicalForm dummy;
     1560      CFList recDiophantine;
     1561      CFList buf2, buf3;
     1562      buf2= factors;
     1563      buf3= p;
     1564      for (CFListIterator iii= buf2; iii.hasItem(); iii++)
     1565        iii.getItem()= mod (iii.getItem(), y);
     1566      for (CFListIterator iii= buf3; iii.hasItem(); iii++)
     1567        iii.getItem()= mod (iii.getItem(), y);
     1568      recDiophantine= diophantine (recResult, buf2, buf3, buf, coeffE);
     1569      CFListIterator iter= recDiophantine;
     1570      for (; j.hasItem(); j++, k++, l++, ii++, iter++)
     1571      {
     1572        k.getItem() += iter.getItem()*power (y, i);
     1573        e -= mulMod (iter.getItem()*power (y, i), l.getItem(), M);
     1574      }
     1575    }
    11591576    if (e.isZero())
    11601577      break;
    … …  
    14101827}
    14111828
     1829// so far only tested for two! factor Hensel lifting
     1830static inline
     1831void
     1832henselStep2 (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
     1833            const CFList& diophant, CFMatrix& M, CFArray& Pi,
     1834            const CFList& products, int j, const CFList& MOD)
     1835{
     1836  CanonicalForm E;
     1837  CanonicalForm xToJ= power (F.mvar(), j);
     1838  Variable x= F.mvar();
     1839
     1840  // compute the error
     1841#ifdef DEBUGOUTPUT
     1842    CanonicalForm test= 1;
     1843    for (int i= 0; i < factors.length(); i++)
     1844    {
     1845      if (i == 0)
     1846        test *= mod (bufFactors [i], power (x, j));
     1847      else
     1848        test *= bufFactors[i];
     1849    }
     1850    test= mod (test, power (x, j));
     1851    test= mod (test, MOD);
     1852    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
     1853    DEBOUTLN (cerr, "test= " << test2);
     1854#endif
     1855
     1856  if (degree (Pi [factors.length() - 2], x) > 0)
     1857    E= F[j] - Pi [factors.length() - 2] [j];
     1858  else
     1859    E= F[j];
     1860
     1861  CFArray buf= CFArray (diophant.length());
     1862
     1863  // actual lifting
     1864  CFList diophantine2= diophantine (diophant, factors, products, MOD, E);
     1865
     1866  int k= 0;
     1867  for (CFListIterator i= diophantine2; k < factors.length(); k++, i++)
     1868  {
     1869    buf[k]= i.getItem();
     1870    bufFactors[k] += xToJ*i.getItem();
     1871  }
     1872
     1873  // update Pi [0]
     1874  int degBuf0= degree (bufFactors[0], x);
     1875  int degBuf1= degree (bufFactors[1], x);
     1876  if (degBuf0 > 0 && degBuf1 > 0)
     1877  {
     1878    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
     1879    if (j + 2 <= M.rows())
     1880      M (j + 2, 1)= mulMod (bufFactors[0] [j + 1], bufFactors[1] [j + 1], MOD);
     1881  }
     1882  CanonicalForm uIZeroJ;
     1883
     1884    if (degBuf0 > 0 && degBuf1 > 0)
     1885      uIZeroJ= mulMod (bufFactors[0] [0], buf[1], MOD) +
     1886               mulMod (bufFactors[1] [0], buf[0], MOD);
     1887    else if (degBuf0 > 0)
     1888      uIZeroJ= mulMod (buf[0], bufFactors[1], MOD);
     1889    else if (degBuf1 > 0)
     1890      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
     1891    else
     1892      uIZeroJ= 0;
     1893    Pi [0] += xToJ*uIZeroJ;
     1894
     1895  CFArray tmp= CFArray (factors.length() - 1);
     1896  for (k= 0; k < factors.length() - 1; k++)
     1897    tmp[k]= 0;
     1898  CFIterator one, two;
     1899  one= bufFactors [0];
     1900  two= bufFactors [1];
     1901  if (degBuf0 > 0 && degBuf1 > 0)
     1902  {
     1903    while (one.exp() > j) one++;
     1904    while (two.exp() > j) two++;
     1905    for (k= 1; k <= (int) ceil (j/2.0); k++)
     1906    { 
     1907      if (k != j - k + 1)
     1908      {
     1909        if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
     1910        {
     1911          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
     1912                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
     1913                    M (j - k + 2, 1);
     1914          one++;
     1915          two++;
     1916        }
     1917        else if (one.exp() == j - k + 1)
     1918        {
     1919          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
     1920                    bufFactors[1] [k], MOD) - M (k + 1, 1);
     1921          one++;
     1922        }
     1923        else if (two.exp() == j - k + 1)
     1924        {
     1925          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
     1926                    two.coeff()), MOD) - M (k + 1, 1);
     1927          two++;
     1928        }
     1929      }
     1930      else
     1931      {
     1932        tmp[0] += M (k + 1, 1);
     1933      }
     1934    }
     1935   
     1936    if (j + 2 <= M.rows())
     1937    {
     1938      if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
     1939        tmp [0] += mulMod ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
     1940                           (bufFactors [1] [j + 1] + bufFactors [1] [0]), MOD)
     1941                   - M(1,1) - M (j + 2,1);
     1942      else if (degBuf0 >= j + 1)
     1943        tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1] [0], MOD);
     1944      else if (degBuf1 >= j + 1)
     1945        tmp[0] += mulMod (bufFactors [0] [0], bufFactors [1] [j + 1], MOD);
     1946    }
     1947  }
     1948  Pi [0] += tmp[0]*xToJ*F.mvar();
     1949
     1950  // update Pi [l]
     1951  int degPi, degBuf;
     1952  for (int l= 1; l < factors.length() - 1; l++)
     1953  {
     1954    degPi= degree (Pi [l - 1], x);
     1955    degBuf= degree (bufFactors[l + 1], x);
     1956    if (degPi > 0 && degBuf > 0)
     1957      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD); 
     1958    if (j == 1)
     1959    {
     1960      if (degPi > 0 && degBuf > 0)
     1961        Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]),
     1962                  (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)-
     1963                  M (1, l + 1));
     1964      else if (degPi > 0)
     1965        Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD));
     1966      else if (degBuf > 0)
     1967        Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD));
     1968    }
     1969    else
     1970    {
     1971      if (degPi > 0 && degBuf > 0)
     1972      {
     1973        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
     1974        uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD);
     1975      }
     1976      else if (degPi > 0)
     1977        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD);
     1978      else if (degBuf > 0)
     1979      {
     1980        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
     1981        uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD);
     1982      }
     1983      Pi[l] += xToJ*uIZeroJ;
     1984    }
     1985    one= bufFactors [l + 1];
     1986    two= Pi [l - 1];
     1987    if (two.exp() == j + 1)
     1988    {
     1989      if (degBuf > 0 && degPi > 0)
     1990      {
     1991          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD);
     1992          two++;
     1993      }
     1994      else if (degPi > 0)
     1995      {     
     1996          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD);
     1997          two++;
     1998      }
     1999    }
     2000    if (degBuf > 0 && degPi > 0)
     2001    {
     2002      for (k= 1; k <= (int) ceil (j/2.0); k++)
     2003      { 
     2004        if (k != j - k + 1)
     2005        {
     2006          if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
     2007          {
     2008            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
     2009                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
     2010                      M (j - k + 2, l + 1);
     2011            one++;
     2012            two++;
     2013          }
     2014          else if (one.exp() == j - k + 1)
     2015          {
     2016            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
     2017                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);     
     2018            one++;
     2019          }
     2020          else if (two.exp() == j - k + 1)
     2021          {
     2022            tmp[l] += mulMod (bufFactors[l + 1] [k],
     2023                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
     2024            two++;
     2025          }
     2026        }
     2027        else
     2028          tmp[l] += M (k + 1, l + 1);
     2029      } 
     2030    } 
     2031    Pi[l] += tmp[l]*xToJ*F.mvar();
     2032  }
     2033
     2034  return;
     2035}
     2036
    14122037CFList
    14132038henselLift23 (const CFList& eval, const CFList& factors, const int* l, CFList&
    … …  
    14162041  CFList buf= factors;
    14172042  int k= 0;
    1418   int liftBound;
    14192043  int liftBoundBivar= l[k];
    14202044  diophant= biDiophantine (eval.getFirst(), buf, liftBoundBivar);
    … …  
    14382062  if (i.hasItem())
    14392063    i++;
    1440   for (i; i.hasItem(); i++, k++)
     2064  for (; i.hasItem(); i++, k++)
    14412065  {
    14422066    Pi [k]= mulMod (Pi [k - 1], i.getItem(), MOD);
    … …  
    15052129  if (i.hasItem())
    15062130    i++;
    1507   for (i; i.hasItem(); i++, k++)
     2131  for (; i.hasItem(); i++, k++)
    15082132  {
    15092133    Pi [k]= mod (Pi [k], xToLOld);
    … …  
    15212145CFList
    15222146henselLift (const CFList& eval, const CFList& factors, const int* l, const int
    1523             lLength)
     2147            lLength, bool sort)
    15242148{
    15252149  CFList diophant;
    15262150  CFList buf= factors;
    15272151  buf.insert (LC (eval.getFirst(), 1));
    1528   sortList (buf, Variable (1));
     2152  if (sort)
     2153    sortList (buf, Variable (1));
    15292154  CFArray Pi;
    15302155  CFMatrix M= CFMatrix (l[1], factors.length());
    … …  
    15532178}
    15542179
     2180// wrt. Variable (1)
     2181CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c)
     2182{
     2183  if (degree (F, 1) <= 0)
     2184   return c;
     2185  else
     2186  {
     2187    CanonicalForm result= swapvar (F, Variable (F.level() + 1), Variable (1));
     2188    result += (swapvar (c, Variable (F.level() + 1), Variable (1))
     2189              - LC (result))*power (result.mvar(), degree (result));
     2190    return swapvar (result, Variable (F.level() + 1), Variable (1));
     2191  }
     2192}
     2193
     2194// so far only for two factor Hensel lifting
     2195CFList
     2196henselLift232 (const CFList& eval, const CFList& factors, const int* l, CFList&
     2197              diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
     2198              const CFList& LCs2)
     2199{
     2200  CFList buf= factors;
     2201  int k= 0;
     2202  int liftBoundBivar= l[k];
     2203  CFList bufbuf= factors;
     2204  Variable v= Variable (2);
     2205
     2206  CFList MOD;
     2207  MOD.append (power (Variable (2), liftBoundBivar));
     2208  CFArray bufFactors= CFArray (factors.length());
     2209  k= 0;
     2210  CFListIterator j= eval;
     2211  j++;
     2212  CFListIterator iter1= LCs1;
     2213  CFListIterator iter2= LCs2;
     2214  iter1++;
     2215  iter2++;
     2216  bufFactors[0]= replaceLC (buf.getFirst(), iter1.getItem());
     2217  bufFactors[1]= replaceLC (buf.getLast(), iter2.getItem());
     2218
     2219  CFListIterator i= buf;
     2220  i++;
     2221  Variable y= j.getItem().mvar();
     2222  if (y.level() != 3)
     2223    y= Variable (3);
     2224
     2225  Pi[0]= mod (Pi[0], power (v, liftBoundBivar));
     2226  M (1, 1)= Pi[0];
     2227  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
     2228    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
     2229                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
     2230  else if (degree (bufFactors[0], y) > 0)
     2231    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
     2232  else if (degree (bufFactors[1], y) > 0)
     2233    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
     2234 
     2235  CFList products;
     2236  for (int i= 0; i < bufFactors.size(); i++)
     2237  {
     2238    if (degree (bufFactors[i], y) > 0)
     2239      products.append (M (1, 1)/bufFactors[i] [0]);
     2240    else
     2241      products.append (M (1, 1)/bufFactors[i]);
     2242  }
     2243 
     2244  for (int d= 1; d < l[1]; d++)
     2245    henselStep2 (j.getItem(), buf, bufFactors, diophant, M, Pi, products, d, MOD);
     2246  CFList result;
     2247  for (k= 0; k < factors.length(); k++)
     2248    result.append (bufFactors[k]);
     2249  return result;
     2250}
     2251
     2252
     2253CFList
     2254henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
     2255            diophant, CFArray& Pi, CFMatrix& M, const int lOld, const int
     2256            lNew, const CFList& LCs1, const CFList& LCs2)
     2257{
     2258  int k= 0;
     2259  CFArray bufFactors= CFArray (factors.length());
     2260  bufFactors[0]= replaceLC (factors.getFirst(), LCs1.getLast());
     2261  bufFactors[1]= replaceLC (factors.getLast(), LCs2.getLast());
     2262  CFList buf= factors;
     2263  Variable y= F.getLast().mvar();
     2264  Variable x= F.getFirst().mvar();
     2265  CanonicalForm xToLOld= power (x, lOld);
     2266  Pi [0]= mod (Pi[0], xToLOld);
     2267  M (1, 1)= Pi [0];
     2268
     2269  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
     2270    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
     2271                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
     2272  else if (degree (bufFactors[0], y) > 0)
     2273    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
     2274  else if (degree (bufFactors[1], y) > 0)
     2275    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
     2276 
     2277  CFList products;
     2278  for (int i= 0; i < bufFactors.size(); i++)
     2279  {
     2280    if (degree (bufFactors[i], y) > 0)
     2281    {
     2282      ASSERT (fdivides (bufFactors[i] [0], M (1, 1)), "expected exact division");
     2283      products.append (M (1, 1)/bufFactors[i] [0]);
     2284    }
     2285    else
     2286    {
     2287      ASSERT (fdivides (bufFactors[i], M (1, 1)), "expected exact division");
     2288      products.append (M (1, 1)/bufFactors[i]);
     2289    }
     2290  }
     2291 
     2292  for (int d= 1; d < lNew; d++)
     2293    henselStep2 (F.getLast(), buf, bufFactors, diophant, M, Pi, products, d, MOD);
     2294  CFList result;
     2295  for (k= 0; k < factors.length(); k++)
     2296    result.append (bufFactors[k]);
     2297  return result;
     2298}
     2299
     2300// so far only for two! factor Hensel lifting
     2301CFList
     2302henselLift2 (const CFList& eval, const CFList& factors, const int* l, const int
     2303             lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
     2304             const CFArray& Pi, const CFList& diophant)
     2305{
     2306  CFList bufDiophant= diophant;
     2307  CFList buf= factors;
     2308  if (sort)
     2309    sortList (buf, Variable (1));
     2310  CFArray bufPi= Pi;
     2311  CFMatrix M= CFMatrix (l[1], factors.length());
     2312  CFList result= henselLift232 (eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2);
     2313  if (eval.length() == 2)
     2314    return result;
     2315  CFList MOD;
     2316  for (int i= 0; i < 2; i++)
     2317    MOD.append (power (Variable (i + 2), l[i]));
     2318  CFListIterator j= eval;
     2319  j++;
     2320  CFList bufEval;
     2321  bufEval.append (j.getItem());
     2322  j++;
     2323  CFListIterator jj= LCs1;
     2324  CFListIterator jjj= LCs2;
     2325  CFList bufLCs1, bufLCs2;
     2326  jj++, jjj++;
     2327  bufLCs1.append (jj.getItem());
     2328  bufLCs2.append (jjj.getItem());
     2329  jj++, jjj++;
     2330 
     2331  for (int i= 2; i <= lLength && j.hasItem(); i++, j++, jj++, jjj++)
     2332  {
     2333    bufEval.append (j.getItem());
     2334    bufLCs1.append (jj.getItem());
     2335    bufLCs2.append (jjj.getItem());
     2336    M= CFMatrix (l[i], factors.length());
     2337    result= henselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M, l[i - 1],
     2338                         l[i], bufLCs1, bufLCs2);
     2339    MOD.append (power (Variable (i + 2), l[i]));
     2340    bufEval.removeFirst();
     2341    bufLCs1.removeFirst();
     2342    bufLCs2.removeFirst();
     2343  }
     2344  return result;
     2345}
     2346
    15552347#endif
    15562348/* HAVE_NTL */

  • factory/facHensel.h

    rce35eb r139b67  
    364364           );
    365365
     366/// two factor Hensel lifting from univariate to bivariate, factors need not to
     367/// be monic
     368void
     369henselLift122 (const CanonicalForm& F,///< [in] a bivariate poly
     370               CFList& factors,       ///< [in, out] a list of univariate polys
     371                                      ///< lifted factors
     372               int l,                 ///< [in] lift bound
     373               CFArray& Pi,           ///< [in, out] stores intermediate results
     374               CFList& diophant,      ///< [in, out] result of diophantine
     375               CFMatrix& M,           ///< [in, out] stores intermediate results
     376               const CFArray& LCs,    ///< [in] leading coefficients
     377               bool sort              ///< [in] if true factors are sorted by
     378                                      ///< their degree
     379              );
     380     
     381
     382/// two factor Hensel lifting from bivariate to multivariate, factors need not
     383/// to be monic
     384///
     385/// @return @a henselLift122 returns a list of lifted factors   
     386CFList
     387henselLift2 (const CFList& eval,   ///< [in] a list of polynomials the last
     388                                   ///< element is a compressed multivariate
     389                                   ///< poly, last but one element equals the
     390                                   ///< last elements modulo its main variable
     391                                   ///< and so on. The first element is a
     392                                   ///< compressed bivariate poly.
     393             const CFList& factors,///< [in] bivariate factors
     394             const int* l,         ///< [in] lift bounds
     395             const int lLength,    ///< [in] length of l
     396             bool sort,            ///< [in] if true factors are sorted by
     397                                   ///< their degree in Variable(1)
     398             const CFList& LCs1,   ///< [in] a list of evaluated LC of first
     399                                   ///< factor
     400             const CFList& LCs2,   ///< [in] a list of evaluated LC of second
     401                                   ///< factor
     402             const CFArray& Pi,    ///< [in] intermediate result
     403             const CFList& diophant///< [in] result of diophantine
     404            );
    366405#endif
    367406/* FAC_HENSEL_H */
Note: See TracChangeset for help on using the changeset viewer.