Changeset 08daea in git for factory/facHensel.cc


Ignore:
Timestamp:
Nov 22, 2010, 11:12:46 AM (13 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
fb8ba2756f35787214f1a0b5702e31e4b853de97
Parents:
cde7084dc141670809afd29cafbd9c571bd9fc3a
Message:
new ezgcd, heuristic which gcd to choose and sparse modular gcd all for finite fields


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

Legend:

Unmodified
Added
Removed

  • factory/facHensel.cc

    rcde708 r08daea  
    2626#include "facHensel.h"
    2727#include "cf_util.h"
     28#include "fac_util.h"
     29#include "cf_algorithm.h"
    2830
    2931#ifdef HAVE_NTL
    … …  
    15531555}
    15541556
     1557void
     1558henselStep122 (const CanonicalForm& F, const CFList& factors,
     1559              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
     1560              CFArray& Pi, int j, const CFArray& LCs)
     1561{
     1562  Variable x= F.mvar();
     1563  CanonicalForm xToJ= power (x, j);
     1564
     1565  CanonicalForm E;
     1566  // compute the error
     1567  if (degree (Pi [factors.length() - 2], x) > 0)
     1568    E= F[j] - Pi [factors.length() - 2] [j];
     1569  else
     1570    E= F[j];
     1571
     1572  CFArray buf= CFArray (diophant.length());
     1573
     1574  int k= 0;
     1575  CanonicalForm remainder;
     1576  // actual lifting
     1577  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
     1578  {
     1579    if (degree (bufFactors[k], x) > 0)
     1580      remainder= modNTL (E, bufFactors[k] [0]);
     1581    else
     1582      remainder= modNTL (E, bufFactors[k]);
     1583    buf[k]= mulNTL (i.getItem(), remainder);
     1584    if (degree (bufFactors[k], x) > 0)
     1585      buf[k]= modNTL (buf[k], bufFactors[k] [0]);
     1586    else
     1587      buf[k]= modNTL (buf[k], bufFactors[k]);
     1588  }
     1589
     1590  for (k= 0; k < factors.length(); k++)
     1591    bufFactors[k] += xToJ*buf[k];
     1592
     1593  // update Pi [0]
     1594  int degBuf0= degree (bufFactors[0], x);
     1595  int degBuf1= degree (bufFactors[1], x);
     1596  if (degBuf0 > 0 && degBuf1 > 0)
     1597  {
     1598    M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
     1599    if (j + 2 <= M.rows())
     1600      M (j + 2, 1)= mulNTL (bufFactors[0] [j + 1], bufFactors[1] [j + 1]);
     1601  }
     1602  CanonicalForm uIZeroJ;
     1603  if (degBuf0 > 0 && degBuf1 > 0)
     1604    uIZeroJ= mulNTL(bufFactors[0][0],buf[1])+mulNTL (bufFactors[1][0], buf[0]);
     1605  else if (degBuf0 > 0)
     1606    uIZeroJ= mulNTL (buf[0], bufFactors[1]);
     1607  else if (degBuf1 > 0)
     1608    uIZeroJ= mulNTL (bufFactors[0], buf [1]);
     1609  else
     1610    uIZeroJ= 0;
     1611  Pi [0] += xToJ*uIZeroJ;
     1612
     1613  CFArray tmp= CFArray (factors.length() - 1);
     1614  for (k= 0; k < factors.length() - 1; k++)
     1615    tmp[k]= 0;
     1616  CFIterator one, two;
     1617  one= bufFactors [0];
     1618  two= bufFactors [1];
     1619  if (degBuf0 > 0 && degBuf1 > 0)
     1620  {
     1621    while (one.exp() > j) one++;
     1622    while (two.exp() > j) two++;
     1623    for (k= 1; k <= (int) ceil (j/2.0); k++)
     1624    {
     1625      if (k != j - k + 1)
     1626      {
     1627        if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
     1628        {
     1629          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()),(bufFactors[1][k] +
     1630                      two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
     1631          one++;
     1632          two++;
     1633        }
     1634        else if (one.exp() == j - k + 1)
     1635        {
     1636          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1] [k]) -
     1637                      M (k + 1, 1);
     1638          one++;
     1639        }
     1640        else if (two.exp() == j - k + 1)
     1641        {
     1642          tmp[0] += mulNTL (bufFactors[0][k],(bufFactors[1][k] + two.coeff())) -
     1643                    M (k + 1, 1);
     1644          two++;
     1645        }
     1646      }
     1647      else
     1648        tmp[0] += M (k + 1, 1);
     1649    }
     1650
     1651    if (j + 2 <= M.rows())
     1652    {
     1653      if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
     1654        tmp [0] += mulNTL ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
     1655                           (bufFactors [1] [j + 1] + bufFactors [1] [0]))
     1656                           - M(1,1) - M (j + 2,1);
     1657      else if (degBuf0 >= j + 1)
     1658        tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1] [0]);
     1659      else if (degBuf1 >= j + 1)
     1660        tmp[0] += mulNTL (bufFactors [0] [0], bufFactors [1] [j + 1]);
     1661    }
     1662  }
     1663  Pi [0] += tmp[0]*xToJ*F.mvar();
     1664
     1665  /*// update Pi [l]
     1666  int degPi, degBuf;
     1667  for (int l= 1; l < factors.length() - 1; l++)
     1668  {
     1669    degPi= degree (Pi [l - 1], x);
     1670    degBuf= degree (bufFactors[l + 1], x);
     1671    if (degPi > 0 && degBuf > 0)
     1672      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
     1673    if (j == 1)
     1674    {
     1675      if (degPi > 0 && degBuf > 0)
     1676        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
     1677                  bufFactors[l + 1] [0] + buf[l + 1]) - M (j + 1, l +1) -
     1678                  M (1, l + 1));
     1679      else if (degPi > 0)
     1680        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1]));
     1681      else if (degBuf > 0)
     1682        Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1]));
     1683    }
     1684    else
     1685    {
     1686      if (degPi > 0 && degBuf > 0)
     1687      {
     1688        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
     1689        uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1]);
     1690      }
     1691      else if (degPi > 0)
     1692        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1]);
     1693      else if (degBuf > 0)
     1694      {
     1695        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
     1696        uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1]);
     1697      }
     1698      Pi[l] += xToJ*uIZeroJ;
     1699    }
     1700    one= bufFactors [l + 1];
     1701    two= Pi [l - 1];
     1702    if (two.exp() == j + 1)
     1703    {
     1704      if (degBuf > 0 && degPi > 0)
     1705      {
     1706          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0]);
     1707          two++;
     1708      }
     1709      else if (degPi > 0)
     1710      {
     1711          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1]);
     1712          two++;
     1713      }
     1714    }
     1715    if (degBuf > 0 && degPi > 0)
     1716    {
     1717      for (k= 1; k <= (int) ceil (j/2.0); k++)
     1718      {
     1719        if (k != j - k + 1)
     1720        {
     1721          if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
     1722          {
     1723            tmp[l] += mulNTL ((bufFactors[l + 1][k]+one.coeff()),(Pi[l-1] [k] +
     1724                        two.coeff())) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
     1725            one++;
     1726            two++;
     1727          }
     1728          else if (one.exp() == j - k + 1)
     1729          {
     1730            tmp[l] += mulNTL ((bufFactors[l+1][k]+one.coeff()),Pi[l - 1] [k]) -
     1731                        M (k + 1, l + 1);
     1732            one++;
     1733          }
     1734          else if (two.exp() == j - k + 1)
     1735          {
     1736            tmp[l] += mulNTL (bufFactors[l+1][k],(Pi[l-1] [k] + two.coeff())) -
     1737                      M (k + 1, l + 1);
     1738            two++;
     1739          }
     1740        }
     1741        else
     1742          tmp[l] += M (k + 1, l + 1);
     1743      }
     1744    }
     1745    Pi[l] += tmp[l]*xToJ*F.mvar();
     1746  }*/
     1747  return;
     1748}
     1749
     1750void
     1751henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
     1752              CFList& diophant, CFMatrix& M, const CFArray& LCs, bool sort)
     1753{
     1754  if (sort)
     1755    sortList (factors, Variable (1));
     1756  Pi= CFArray (factors.length() - 2);
     1757  CFList bufFactors2= factors;
     1758  bufFactors2.removeFirst();
     1759  diophant.removeFirst();
     1760  CFListIterator iter= diophant;
     1761  CanonicalForm s,t;
     1762  extgcd (bufFactors2.getFirst(), bufFactors2.getLast(), s, t);
     1763  diophant= CFList();
     1764  diophant.append (t);
     1765  diophant.append (s);
     1766  DEBOUTLN (cerr, "diophant= " << diophant);
     1767
     1768  CFArray bufFactors= CFArray (bufFactors2.length());
     1769  int i= 0;
     1770  for (CFListIterator k= bufFactors2; k.hasItem(); i++, k++)
     1771    bufFactors[i]= replaceLc (k.getItem(), LCs [i]);
     1772
     1773  Variable x= F.mvar();
     1774  if (degree (bufFactors[0], x) > 0 && degree (bufFactors [1], x) > 0)
     1775  {
     1776    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1] [0]);
     1777    Pi [0]= M (1, 1) + (mulNTL (bufFactors [0] [1], bufFactors[1] [0]) +
     1778                        mulNTL (bufFactors [0] [0], bufFactors [1] [1]))*x;
     1779  }
     1780  else if (degree (bufFactors[0], x) > 0)
     1781  {
     1782    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1]);
     1783    Pi [0]= M (1, 1) +
     1784            mulNTL (bufFactors [0] [1], bufFactors[1])*x;
     1785  }
     1786  else if (degree (bufFactors[1], x) > 0)
     1787  {
     1788    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1] [0]);
     1789    Pi [0]= M (1, 1) +
     1790            mulNTL (bufFactors [0], bufFactors[1] [1])*x;
     1791  }
     1792  else
     1793  {
     1794    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1]);
     1795    Pi [0]= M (1, 1);
     1796  }
     1797
     1798  for (i= 1; i < l; i++)
     1799    henselStep122 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
     1800
     1801  factors= CFList();
     1802  for (i= 0; i < bufFactors.size(); i++)
     1803    factors.append (bufFactors[i]);
     1804  return;
     1805}
     1806
     1807/// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$
     1808/// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$
     1809static inline
     1810CFList
     1811diophantine (const CFList& recResult, const CFList& factors,
     1812             const CFList& products, const CFList& M, const CanonicalForm& E)
     1813{
     1814  if (M.isEmpty())
     1815  {
     1816    CFList result;
     1817    CFListIterator j= factors;
     1818    CanonicalForm buf;
     1819    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
     1820    {
     1821      ASSERT (E.isUnivariate() || E.inCoeffDomain(),
     1822              "constant or univariate poly expected");
     1823      ASSERT (i.getItem().isUnivariate() || i.getItem().inCoeffDomain(),
     1824              "constant or univariate poly expected");
     1825      ASSERT (j.getItem().isUnivariate() || j.getItem().inCoeffDomain(),
     1826              "constant or univariate poly expected");
     1827      buf= mulNTL (E, i.getItem());
     1828      result.append (modNTL (buf, j.getItem()));
     1829    }
     1830    return result;
     1831  }
     1832  Variable y= M.getLast().mvar();
     1833  CFList bufFactors= factors;
     1834  for (CFListIterator i= bufFactors; i.hasItem(); i++)
     1835    i.getItem()= mod (i.getItem(), y);
     1836  CFList bufProducts= products;
     1837  for (CFListIterator i= bufProducts; i.hasItem(); i++)
     1838    i.getItem()= mod (i.getItem(), y);
     1839  CFList buf= M;
     1840  buf.removeLast();
     1841  CanonicalForm bufE= mod (E, y);
     1842  CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
     1843                                      bufE);
     1844
     1845  CanonicalForm e= E;
     1846  CFListIterator j= products;
     1847  for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++)
     1848    e -= mulMod (i.getItem(), j.getItem(), M);
     1849
     1850  CFList result= recDiophantine;
     1851  int d= degree (M.getLast());
     1852  CanonicalForm coeffE;
     1853  for (int i= 1; i < d; i++)
     1854  {
     1855    if (degree (e, y) > 0)
     1856      coeffE= e[i];
     1857    else
     1858      coeffE= 0;
     1859    if (!coeffE.isZero())
     1860    {
     1861      CFListIterator k= result;
     1862      recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
     1863                                   coeffE);
     1864      CFListIterator l= products;
     1865      for (j= recDiophantine; j.hasItem(); j++, k++, l++)
     1866      {
     1867        k.getItem() += j.getItem()*power (y, i);
     1868        e -= mulMod (j.getItem()*power (y, i), l.getItem(), M);
     1869      }
     1870    }
     1871    if (e.isZero())
     1872      break;
     1873  }
     1874  return result;
     1875}
     1876
     1877// non monic case
     1878static inline
     1879CFList
     1880multiRecDiophantine2 (const CanonicalForm& F, const CFList& factors,
     1881                     const CFList& recResult, const CFList& M, const int d)
     1882{
     1883  Variable y= F.mvar();
     1884  CFList result;
     1885  CFListIterator i;
     1886  CanonicalForm e= 1;
     1887  CFListIterator j= factors;
     1888  CFList p;
     1889  CFArray bufFactors= CFArray (factors.length());
     1890  CanonicalForm yToD= power (y, d);
     1891  int k= 0;
     1892  for (CFListIterator i= factors; i.hasItem(); i++, k++)
     1893    bufFactors [k]= i.getItem();
     1894  CanonicalForm b;
     1895  CFList buf= M;
     1896  buf.removeLast();
     1897  buf.append (yToD);
     1898  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
     1899  {
     1900    b= 1;
     1901    if (fdivides (bufFactors[k], F))
     1902      b= F/bufFactors[k];
     1903    else
     1904    {
     1905      for (int l= 0; l < factors.length(); l++)
     1906      {
     1907        if (l == k)
     1908          continue;
     1909        else
     1910        {
     1911          b= mulMod (b, bufFactors[l], buf);
     1912        }
     1913      }
     1914    }
     1915    p.append (b);
     1916  }
     1917  j= p;
     1918  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
     1919    e -= mulMod (i.getItem(), j.getItem(), M);
     1920  if (e.isZero())
     1921    return recResult;
     1922  CanonicalForm coeffE;
     1923  result= recResult;
     1924  CanonicalForm g;
     1925  buf.removeLast();
     1926  for (int i= 1; i < d; i++)
     1927  {
     1928    if (degree (e, y) > 0)
     1929      coeffE= e[i];
     1930    else
     1931      coeffE= 0;
     1932    if (!coeffE.isZero())
     1933    {
     1934      CFListIterator k= result;
     1935      CFListIterator l= p;
     1936      j= recResult;
     1937      int ii= 0;
     1938      CanonicalForm dummy;
     1939      CFList recDiophantine;
     1940      CFList buf2, buf3;
     1941      buf2= factors;
     1942      buf3= p;
     1943      for (CFListIterator iii= buf2; iii.hasItem(); iii++)
     1944        iii.getItem()= mod (iii.getItem(), y);
     1945      for (CFListIterator iii= buf3; iii.hasItem(); iii++)
     1946        iii.getItem()= mod (iii.getItem(), y);
     1947      recDiophantine= diophantine (recResult, buf2, buf3, buf, coeffE);
     1948      CFListIterator iter= recDiophantine;
     1949      for (; j.hasItem(); j++, k++, l++, ii++, iter++)
     1950      {
     1951        k.getItem() += iter.getItem()*power (y, i);
     1952        e -= mulMod (iter.getItem()*power (y, i), l.getItem(), M);
     1953      }
     1954    }
     1955    if (e.isZero())
     1956      break;
     1957  }
     1958
     1959#ifdef DEBUGOUTPUT
     1960    CanonicalForm test= 0;
     1961    j= p;
     1962    for (CFListIterator i= result; i.hasItem(); i++, j++)
     1963      test += mod (i.getItem()*j.getItem(), power (y, d));
     1964    DEBOUTLN (cerr, "test= " << test);
     1965#endif
     1966  return result;
     1967}
     1968
     1969// so far only tested for two! factor Hensel lifting
     1970static inline
     1971void
     1972henselStep2 (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
     1973            const CFList& diophant, CFMatrix& M, CFArray& Pi,
     1974            const CFList& products, int j, const CFList& MOD)
     1975{
     1976  CanonicalForm E;
     1977  CanonicalForm xToJ= power (F.mvar(), j);
     1978  Variable x= F.mvar();
     1979
     1980  // compute the error
     1981#ifdef DEBUGOUTPUT
     1982    CanonicalForm test= 1;
     1983    for (int i= 0; i < factors.length(); i++)
     1984    {
     1985      if (i == 0)
     1986        test *= mod (bufFactors [i], power (x, j));
     1987      else
     1988        test *= bufFactors[i];
     1989    }
     1990    test= mod (test, power (x, j));
     1991    test= mod (test, MOD);
     1992    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
     1993    DEBOUTLN (cerr, "test= " << test2);
     1994#endif
     1995
     1996  if (degree (Pi [factors.length() - 2], x) > 0)
     1997    E= F[j] - Pi [factors.length() - 2] [j];
     1998  else
     1999    E= F[j];
     2000
     2001  CFArray buf= CFArray (diophant.length());
     2002
     2003  // actual lifting
     2004  CFList diophantine2= diophantine (diophant, factors, products, MOD, E);
     2005
     2006  int k= 0;
     2007  for (CFListIterator i= diophantine2; k < factors.length(); k++, i++)
     2008  {
     2009    buf[k]= i.getItem();
     2010    bufFactors[k] += xToJ*i.getItem();
     2011  }
     2012
     2013  // update Pi [0]
     2014  int degBuf0= degree (bufFactors[0], x);
     2015  int degBuf1= degree (bufFactors[1], x);
     2016  if (degBuf0 > 0 && degBuf1 > 0)
     2017  {
     2018    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
     2019    if (j + 2 <= M.rows())
     2020      M (j + 2, 1)= mulMod (bufFactors[0] [j + 1], bufFactors[1] [j + 1], MOD);
     2021  }
     2022  CanonicalForm uIZeroJ;
     2023
     2024    if (degBuf0 > 0 && degBuf1 > 0)
     2025      uIZeroJ= mulMod (bufFactors[0] [0], buf[1], MOD) +
     2026               mulMod (bufFactors[1] [0], buf[0], MOD);
     2027    else if (degBuf0 > 0)
     2028      uIZeroJ= mulMod (buf[0], bufFactors[1], MOD);
     2029    else if (degBuf1 > 0)
     2030      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
     2031    else
     2032      uIZeroJ= 0;
     2033    Pi [0] += xToJ*uIZeroJ;
     2034
     2035  CFArray tmp= CFArray (factors.length() - 1);
     2036  for (k= 0; k < factors.length() - 1; k++)
     2037    tmp[k]= 0;
     2038  CFIterator one, two;
     2039  one= bufFactors [0];
     2040  two= bufFactors [1];
     2041  if (degBuf0 > 0 && degBuf1 > 0)
     2042  {
     2043    while (one.exp() > j) one++;
     2044    while (two.exp() > j) two++;
     2045    for (k= 1; k <= (int) ceil (j/2.0); k++)
     2046    {
     2047      if (k != j - k + 1)
     2048      {
     2049        if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
     2050        {
     2051          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
     2052                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
     2053                    M (j - k + 2, 1);
     2054          one++;
     2055          two++;
     2056        }
     2057        else if (one.exp() == j - k + 1)
     2058        {
     2059          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
     2060                    bufFactors[1] [k], MOD) - M (k + 1, 1);
     2061          one++;
     2062        }
     2063        else if (two.exp() == j - k + 1)
     2064        {
     2065          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
     2066                    two.coeff()), MOD) - M (k + 1, 1);
     2067          two++;
     2068        }
     2069      }
     2070      else
     2071      {
     2072        tmp[0] += M (k + 1, 1);
     2073      }
     2074    }
     2075
     2076    if (j + 2 <= M.rows())
     2077    {
     2078      if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
     2079        tmp [0] += mulMod ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
     2080                           (bufFactors [1] [j + 1] + bufFactors [1] [0]), MOD)
     2081                   - M(1,1) - M (j + 2,1);
     2082      else if (degBuf0 >= j + 1)
     2083        tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1] [0], MOD);
     2084      else if (degBuf1 >= j + 1)
     2085        tmp[0] += mulMod (bufFactors [0] [0], bufFactors [1] [j + 1], MOD);
     2086    }
     2087  }
     2088  Pi [0] += tmp[0]*xToJ*F.mvar();
     2089
     2090  // update Pi [l]
     2091  int degPi, degBuf;
     2092  for (int l= 1; l < factors.length() - 1; l++)
     2093  {
     2094    degPi= degree (Pi [l - 1], x);
     2095    degBuf= degree (bufFactors[l + 1], x);
     2096    if (degPi > 0 && degBuf > 0)
     2097      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
     2098    if (j == 1)
     2099    {
     2100      if (degPi > 0 && degBuf > 0)
     2101        Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]),
     2102                  (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)-
     2103                  M (1, l + 1));
     2104      else if (degPi > 0)
     2105        Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD));
     2106      else if (degBuf > 0)
     2107        Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD));
     2108    }
     2109    else
     2110    {
     2111      if (degPi > 0 && degBuf > 0)
     2112      {
     2113        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
     2114        uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD);
     2115      }
     2116      else if (degPi > 0)
     2117        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD);
     2118      else if (degBuf > 0)
     2119      {
     2120        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
     2121        uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD);
     2122      }
     2123      Pi[l] += xToJ*uIZeroJ;
     2124    }
     2125    one= bufFactors [l + 1];
     2126    two= Pi [l - 1];
     2127    if (two.exp() == j + 1)
     2128    {
     2129      if (degBuf > 0 && degPi > 0)
     2130      {
     2131          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD);
     2132          two++;
     2133      }
     2134      else if (degPi > 0)
     2135      {
     2136          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD);
     2137          two++;
     2138      }
     2139    }
     2140    if (degBuf > 0 && degPi > 0)
     2141    {
     2142      for (k= 1; k <= (int) ceil (j/2.0); k++)
     2143      {
     2144        if (k != j - k + 1)
     2145        {
     2146          if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
     2147          {
     2148            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
     2149                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
     2150                      M (j - k + 2, l + 1);
     2151            one++;
     2152            two++;
     2153          }
     2154          else if (one.exp() == j - k + 1)
     2155          {
     2156            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
     2157                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
     2158            one++;
     2159          }
     2160          else if (two.exp() == j - k + 1)
     2161          {
     2162            tmp[l] += mulMod (bufFactors[l + 1] [k],
     2163                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
     2164            two++;
     2165          }
     2166        }
     2167        else
     2168          tmp[l] += M (k + 1, l + 1);
     2169      }
     2170    }
     2171    Pi[l] += tmp[l]*xToJ*F.mvar();
     2172  }
     2173
     2174  return;
     2175}
     2176
     2177// wrt. Variable (1)
     2178CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c)
     2179{
     2180  if (degree (F, 1) <= 0)
     2181   return c;
     2182  else
     2183  {
     2184    CanonicalForm result= swapvar (F, Variable (F.level() + 1), Variable (1));
     2185    result += (swapvar (c, Variable (F.level() + 1), Variable (1))
     2186              - LC (result))*power (result.mvar(), degree (result));
     2187    return swapvar (result, Variable (F.level() + 1), Variable (1));
     2188  }
     2189}
     2190
     2191// so far only for two factor Hensel lifting
     2192CFList
     2193henselLift232 (const CFList& eval, const CFList& factors, int* l, CFList&
     2194              diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
     2195              const CFList& LCs2)
     2196{
     2197  CFList buf= factors;
     2198  int k= 0;
     2199  int liftBoundBivar= l[k];
     2200  CFList bufbuf= factors;
     2201  Variable v= Variable (2);
     2202
     2203  CFList MOD;
     2204  MOD.append (power (Variable (2), liftBoundBivar));
     2205  CFArray bufFactors= CFArray (factors.length());
     2206  k= 0;
     2207  CFListIterator j= eval;
     2208  j++;
     2209  CFListIterator iter1= LCs1;
     2210  CFListIterator iter2= LCs2;
     2211  iter1++;
     2212  iter2++;
     2213  bufFactors[0]= replaceLC (buf.getFirst(), iter1.getItem());
     2214  bufFactors[1]= replaceLC (buf.getLast(), iter2.getItem());
     2215
     2216  CFListIterator i= buf;
     2217  i++;
     2218  Variable y= j.getItem().mvar();
     2219  if (y.level() != 3)
     2220    y= Variable (3);
     2221
     2222  Pi[0]= mod (Pi[0], power (v, liftBoundBivar));
     2223  M (1, 1)= Pi[0];
     2224  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
     2225    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
     2226                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
     2227  else if (degree (bufFactors[0], y) > 0)
     2228    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
     2229  else if (degree (bufFactors[1], y) > 0)
     2230    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
     2231
     2232  CFList products;
     2233  for (int i= 0; i < bufFactors.size(); i++)
     2234  {
     2235    if (degree (bufFactors[i], y) > 0)
     2236      products.append (M (1, 1)/bufFactors[i] [0]);
     2237    else
     2238      products.append (M (1, 1)/bufFactors[i]);
     2239  }
     2240
     2241  for (int d= 1; d < l[1]; d++)
     2242    henselStep2 (j.getItem(), buf, bufFactors, diophant, M, Pi, products, d, MOD);
     2243  CFList result;
     2244  for (k= 0; k < factors.length(); k++)
     2245    result.append (bufFactors[k]);
     2246  return result;
     2247}
     2248
     2249
     2250CFList
     2251henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
     2252            diophant, CFArray& Pi, CFMatrix& M, const int lOld, int&
     2253            lNew, const CFList& LCs1, const CFList& LCs2)
     2254{
     2255  int k= 0;
     2256  CFArray bufFactors= CFArray (factors.length());
     2257  bufFactors[0]= replaceLC (factors.getFirst(), LCs1.getLast());
     2258  bufFactors[1]= replaceLC (factors.getLast(), LCs2.getLast());
     2259  CFList buf= factors;
     2260  Variable y= F.getLast().mvar();
     2261  Variable x= F.getFirst().mvar();
     2262  CanonicalForm xToLOld= power (x, lOld);
     2263  Pi [0]= mod (Pi[0], xToLOld);
     2264  M (1, 1)= Pi [0];
     2265
     2266  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
     2267    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
     2268                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
     2269  else if (degree (bufFactors[0], y) > 0)
     2270    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
     2271  else if (degree (bufFactors[1], y) > 0)
     2272    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
     2273
     2274  CFList products;
     2275  for (int i= 0; i < bufFactors.size(); i++)
     2276  {
     2277    if (degree (bufFactors[i], y) > 0)
     2278    {
     2279      ASSERT (fdivides (bufFactors[i][0], M (1, 1)), "expected exact division");
     2280      products.append (M (1, 1)/bufFactors[i] [0]);
     2281    }
     2282    else
     2283    {
     2284      ASSERT (fdivides (bufFactors[i], M (1, 1)), "expected exact division");
     2285      products.append (M (1, 1)/bufFactors[i]);
     2286    }
     2287  }
     2288
     2289  for (int d= 1; d < lNew; d++)
     2290    henselStep2 (F.getLast(),buf,bufFactors, diophant, M, Pi, products, d, MOD);
     2291
     2292  CFList result;
     2293  for (k= 0; k < factors.length(); k++)
     2294    result.append (bufFactors[k]);
     2295  return result;
     2296}
     2297
     2298// so far only for two! factor Hensel lifting
     2299CFList
     2300henselLift2 (const CFList& eval, const CFList& factors, int* l, const int
     2301             lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
     2302             const CFArray& Pi, const CFList& diophant)
     2303{
     2304  CFList bufDiophant= diophant;
     2305  CFList buf= factors;
     2306  if (sort)
     2307    sortList (buf, Variable (1));
     2308  CFArray bufPi= Pi;
     2309  CFMatrix M= CFMatrix (l[1], factors.length());
     2310  CFList result= henselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2);
     2311  if (eval.length() == 2)
     2312    return result;
     2313  CFList MOD;
     2314  for (int i= 0; i < 2; i++)
     2315    MOD.append (power (Variable (i + 2), l[i]));
     2316  CFListIterator j= eval;
     2317  j++;
     2318  CFList bufEval;
     2319  bufEval.append (j.getItem());
     2320  j++;
     2321  CFListIterator jj= LCs1;
     2322  CFListIterator jjj= LCs2;
     2323  CFList bufLCs1, bufLCs2;
     2324  jj++, jjj++;
     2325  bufLCs1.append (jj.getItem());
     2326  bufLCs2.append (jjj.getItem());
     2327  jj++, jjj++;
     2328
     2329  for (int i= 2; i <= lLength && j.hasItem(); i++, j++, jj++, jjj++)
     2330  {
     2331    bufEval.append (j.getItem());
     2332    bufLCs1.append (jj.getItem());
     2333    bufLCs2.append (jjj.getItem());
     2334    M= CFMatrix (l[i], factors.length());
     2335    result= henselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M, l[i - 1],
     2336                         l[i], bufLCs1, bufLCs2);
     2337    MOD.append (power (Variable (i + 2), l[i]));
     2338    bufEval.removeFirst();
     2339    bufLCs1.removeFirst();
     2340    bufLCs2.removeFirst();
     2341  }
     2342  return result;
     2343}
     2344
    15552345#endif
    15562346/* HAVE_NTL */
Note: See TracChangeset for help on using the changeset viewer.