Changeset c4f4fd in git for factory/facHensel.cc


Ignore:
Timestamp:
Sep 16, 2010, 5:34:13 PM (14 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
8e7db4cc8c23b1a45ca71f3217a02819e04d9c84
Parents:
7b8a4166957e36a8fb67484b68a9449846e1995d
Message:
going back to r13182 (facstd.tst fails)

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

Legend:

Unmodified
Added
Removed

  • factory/facHensel.cc

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