Changeset 0e6668 in git


Ignore:
Timestamp:
Jan 6, 2012, 6:21:28 PM (11 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'spielwiese', 'f6c3dc58b0df4bd712574325fe76d0626174ad97')
Children:
84dcd6b39d6c53ee61bc78ee8970dd3f30f0e458
Parents:
5302958e38bf485567be744d126aafa77e73b30e
git-author:
Martin Lee <martinlee84@web.de>2012-01-06 18:21:28+01:00
git-committer:
Martin Lee <martinlee84@web.de>2012-04-04 14:42:25+02:00
Message:
chg: more changes in order to replace old factorization over integers
Location:
factory
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • factory/facBivar.cc

    r530295 r0e6668  
    527527  }
    528528
     529  if (!extension)
     530  {
     531    for (CFListIterator i= uniFactors; i.hasItem(); i++)
     532      i.getItem() /= lc (i.getItem());
     533  }
     534
    529535  A= A (y + evaluation, y);
    530536
     
    534540  CFList earlyFactors;
    535541  TIMING_START (fac_hensel_lift);
     542  //out_cf ("A before= ",A, "\n");
    536543  uniFactors= henselLiftAndEarly0
    537544             (A, earlySuccess, earlyFactors, degs, liftBound,
     
    540547  DEBOUTLN (cerr, "lifted factors= " << uniFactors);
    541548
     549  //printf ("earlyFactors.length()= %d\n", earlyFactors.length());
     550  //printf ("liftBound after= %d\n", liftBound);
     551  //printf ("earlySuccess= %d\n", earlySuccess);
    542552  CanonicalForm MODl= power (y, liftBound);
    543553
     554  CanonicalForm test= prod (uniFactors);
     555  test *= LC (A,1);
     556  test= mod (test, MODl);
     557  //printf ("test == A %d\n", test == A);
     558  //out_cf ("test= ", test, "\n");
     559  //out_cf ("A= ", A, "\n");
    544560  factors= factorRecombination (uniFactors, A, MODl, degs, 1,
    545561                                uniFactors.length()/2);
  • factory/facFactorize.cc

    r530295 r0e6668  
    603603  if (F.isUnivariate())
    604604  {
    605     CFList result= conv (factorize (F, v));
     605    CFList result;
     606    if (v.level() != 1)
     607      result= conv (factorize (F, v));
     608    else
     609      result= conv (factorize (F,true));
    606610    if (result.getFirst().inCoeffDomain())
    607611      result.removeFirst();
     
    654658  if (A.isUnivariate ())
    655659  {
    656     factors= conv (factorize (A, v));
     660    if (v.level() != 1)
     661      factors= conv (factorize (A, v));
     662    else
     663      factors= conv (factorize (A,true));
    657664    append (factors, contentAFactors);
    658665    decompress (factors, N);
  • factory/facHensel.cc

    r530295 r0e6668  
    564564    {
    565565      CanonicalForm result= mulFLINTQa (F, G, alpha);
     566      CanonicalForm test= F*G;
     567      if (test != result)
     568        printf ("FAILLLLLLLL\n");
    566569      return result;
    567570    }
    568571    else if (!F.inCoeffDomain() && !G.inCoeffDomain())
     572    {
     573      CanonicalForm result= mulFLINTQ (F,G);
     574      CanonicalForm test= F*G;
     575      if (test != result)
     576        printf ("FAILLLLLLLL2\n");
    569577      return mulFLINTQ (F, G);
     578    }
    570579#endif
    571580    return F*G;
     
    620629  {
    621630#ifdef HAVE_FLINT
    622     return modFLINTQ (F, G);
     631    Variable alpha;
     632    if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
     633    {
     634      CanonicalForm result= modFLINTQ (F,G);
     635      CanonicalForm test= mod (F, G);
     636      if (test != result)
     637        printf ("FAILINMOD\n");
     638      return result;
     639      //return modFLINTQ (F, G);
     640    }
     641    else
     642     return mod (F, G);
    623643#else
    624644    return mod (F, G);
     
    675695  {
    676696#ifdef HAVE_FLINT
    677     return divFLINTQ (F,G);
     697    Variable alpha;
     698    if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
     699    {
     700      CanonicalForm result= divFLINTQ (F,G);
     701      CanonicalForm test= div (F, G);
     702      if (test != result)
     703        printf ("FAILINDIV\n");
     704      return result;
     705      //return divFLINTQ (F,G);
     706    }
     707    else
     708      return div (F, G);
    678709#else
    679710    return div (F, G);
     
    14121443}
    14131444
     1445#ifdef HAVE_FLINT
     1446CanonicalForm
     1447mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
     1448                CanonicalForm& M);
     1449#endif
     1450
    14141451CanonicalForm
    14151452mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
     
    14491486  }
    14501487  else
     1488#ifdef HAVE_FLINT
     1489    return mulMod2FLINTFp (A, B, M);
     1490#else
    14511491    return mulMod2NTLFp (A, B, M);
    1452 }
     1492#endif
     1493}
     1494
     1495#ifdef HAVE_FLINT
     1496void
     1497kronSubRecipro (nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm& A, int d)
     1498{
     1499  int degAy= degree (A);
     1500  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
     1501  nmod_poly_init2_preinv (subA1, getCharacteristic(), ninv, d*(degAy + 2));
     1502  nmod_poly_init2_preinv (subA2, getCharacteristic(), ninv, d*(degAy + 2));
     1503
     1504  nmod_poly_t buf;
     1505
     1506  for (CFIterator i= A; i.hasTerms(); i++)
     1507  {
     1508    convertFacCF2nmod_poly_t (buf, i.coeff());
     1509
     1510    int k= i.exp()*d;
     1511    int kk= (degAy - i.exp())*d;
     1512    int bufRepLength= (int) nmod_poly_length (buf);
     1513    for (int j= 0; j < bufRepLength; j++)
     1514    {
     1515      nmod_poly_set_coeff_ui (subA1, j + k, n_addmod (nmod_poly_get_coeff_ui (subA1, j+k), nmod_poly_get_coeff_ui (buf, j), getCharacteristic())); //TODO muß man subA1 coeffs von subA1 vorher auf null setzen??
     1516      nmod_poly_set_coeff_ui (subA2, j + kk, n_addmod (nmod_poly_get_coeff_ui (subA2, j + kk),nmod_poly_get_coeff_ui (buf, j), getCharacteristic()));
     1517    }
     1518    nmod_poly_clear (buf);
     1519  }
     1520  _nmod_poly_normalise (subA1);
     1521  _nmod_poly_normalise (subA2);
     1522}
     1523
     1524void
     1525kronSubRecipro (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A, int d)
     1526{
     1527  int degAy= degree (A);
     1528  fmpz_poly_init2 (subA1, d*(degAy + 2));
     1529  fmpz_poly_init2 (subA2, d*(degAy + 2));
     1530
     1531  fmpz_poly_t buf;
     1532  fmpz_t coeff1, coeff2;
     1533
     1534  for (CFIterator i= A; i.hasTerms(); i++)
     1535  {
     1536    convertFacCF2Fmpz_poly_t (buf, i.coeff());
     1537
     1538    int k= i.exp()*d;
     1539    int kk= (degAy - i.exp())*d;
     1540    int bufRepLength= (int) fmpz_poly_length (buf);
     1541    for (int j= 0; j < bufRepLength; j++)
     1542    {
     1543      fmpz_poly_get_coeff_fmpz (coeff1, subA1, j+k); //TODO muß man subA1 coeffs von subA1 vorher auf null setzen??
     1544                                                     // vielleicht ist es schneller fmpz_poly_get_ptr zu nehmen
     1545      fmpz_poly_get_coeff_fmpz (coeff2, buf, j);
     1546      fmpz_add (coeff1, coeff1, coeff2);
     1547      fmpz_poly_set_coeff_fmpz (subA1, j + k, coeff1);
     1548      fmpz_poly_get_coeff_fmpz (coeff1, subA2, j + kk);
     1549      fmpz_add (coeff1, coeff1, coeff2);
     1550      fmpz_poly_set_coeff_fmpz (subA2, j + kk, coeff1);
     1551    }
     1552    fmpz_poly_clear (buf);
     1553  }
     1554  fmpz_clear (coeff1);
     1555  fmpz_clear (coeff2);
     1556  _fmpz_poly_normalise (subA1);
     1557  _fmpz_poly_normalise (subA2);
     1558}
     1559
     1560CanonicalForm reverseSubstQ (const fmpz_poly_t F, int d)
     1561{
     1562  Variable y= Variable (2);
     1563  Variable x= Variable (1);
     1564
     1565  fmpz_poly_t f;
     1566  fmpz_poly_init (f);
     1567  fmpz_poly_set (f, F);
     1568
     1569  fmpz_poly_t buf;
     1570  CanonicalForm result= 0;
     1571  int i= 0;
     1572  int degf= fmpz_poly_degree(f);
     1573  int k= 0;
     1574  int degfSubK;
     1575  int repLength;
     1576  fmpz_t coeff;
     1577  while (degf >= k)
     1578  {
     1579    degfSubK= degf - k;
     1580    if (degfSubK >= d)
     1581      repLength= d;
     1582    else
     1583      repLength= degfSubK + 1;
     1584
     1585    fmpz_poly_init2 (buf, repLength);
     1586    fmpz_init (coeff);
     1587    for (int j= 0; j < repLength; j++)
     1588    {
     1589      fmpz_poly_get_coeff_fmpz (coeff, f, j + k);
     1590      fmpz_poly_set_coeff_fmpz (buf, j, coeff);
     1591    }
     1592    _fmpz_poly_normalise (buf);
     1593
     1594    result += convertFmpz_poly_t2FacCF (buf, x)*power (y, i);
     1595    i++;
     1596    k= d*i;
     1597    fmpz_poly_clear (buf);
     1598    fmpz_clear (coeff);
     1599  }
     1600  fmpz_poly_clear (f);
     1601
     1602  return result;
     1603}
     1604
     1605CanonicalForm
     1606reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
     1607{
     1608  Variable y= Variable (2);
     1609  Variable x= Variable (1);
     1610
     1611  nmod_poly_t f, g;
     1612  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
     1613  nmod_poly_init_preinv (f, getCharacteristic(), ninv);
     1614  nmod_poly_init_preinv (g, getCharacteristic(), ninv);
     1615  nmod_poly_set (f, F);
     1616  nmod_poly_set (g, G);
     1617  int degf= nmod_poly_degree(f);
     1618  int degg= nmod_poly_degree(g);
     1619
     1620
     1621  nmod_poly_t buf1,buf2, buf3;
     1622  /*nmod_poly_init_preinv (buf1, getCharacteristic(), ninv);
     1623  nmod_poly_init_preinv (buf2, getCharacteristic(), ninv);
     1624  nmod_poly_init_preinv (buf3, getCharacteristic(), ninv);*/
     1625
     1626  if (nmod_poly_length (f) < (long) d*(k+1)) //zero padding
     1627    nmod_poly_fit_length (f,(long)d*(k+1));
     1628
     1629  CanonicalForm result= 0;
     1630  int i= 0;
     1631  int lf= 0;
     1632  int lg= d*k;
     1633  int degfSubLf= degf;
     1634  int deggSubLg= degg-lg;
     1635  int repLengthBuf2;
     1636  int repLengthBuf1;
     1637  //zz_p zzpZero= zz_p();
     1638  while (degf >= lf || lg >= 0)
     1639  {
     1640    if (degfSubLf >= d)
     1641      repLengthBuf1= d;
     1642    else if (degfSubLf < 0)
     1643      repLengthBuf1= 0;
     1644    else
     1645      repLengthBuf1= degfSubLf + 1;
     1646    //cout << "repLengthBuf1= " << repLengthBuf1 << "\n";
     1647    //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
     1648    nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
     1649    //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
     1650
     1651    for (int ind= 0; ind < repLengthBuf1; ind++)
     1652      nmod_poly_set_coeff_ui (buf1, ind, nmod_poly_get_coeff_ui (f, ind+lf));
     1653    _nmod_poly_normalise (buf1);
     1654    //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
     1655
     1656    repLengthBuf1= nmod_poly_length (buf1);
     1657
     1658    if (deggSubLg >= d - 1)
     1659      repLengthBuf2= d - 1;
     1660    else if (deggSubLg < 0)
     1661      repLengthBuf2= 0;
     1662    else
     1663      repLengthBuf2= deggSubLg + 1;
     1664
     1665    //cout << "repLengthBuf2= " << repLengthBuf2 << "\n";
     1666    //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
     1667    nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
     1668    //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
     1669    for (int ind= 0; ind < repLengthBuf2; ind++)
     1670      nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
     1671
     1672    _nmod_poly_normalise (buf2);
     1673    //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
     1674    repLengthBuf2= nmod_poly_length (buf2);
     1675
     1676    nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
     1677    for (int ind= 0; ind < repLengthBuf1; ind++)
     1678      nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
     1679    for (int ind= repLengthBuf1; ind < d; ind++)
     1680      nmod_poly_set_coeff_ui (buf3, ind, 0);
     1681    for (int ind= 0; ind < repLengthBuf2; ind++)
     1682      nmod_poly_set_coeff_ui (buf3, ind + d, nmod_poly_get_coeff_ui (buf2, ind));
     1683    _nmod_poly_normalise (buf3);
     1684
     1685    result += convertnmod_poly_t2FacCF (buf3, x)*power (y, i);
     1686    i++;
     1687
     1688
     1689    lf= i*d;
     1690    degfSubLf= degf - lf;
     1691
     1692    lg= d*(k-i);
     1693    deggSubLg= degg - lg;
     1694
     1695    if (lg >= 0 && deggSubLg > 0)
     1696    {
     1697      if (repLengthBuf2 > degfSubLf + 1)
     1698        degfSubLf= repLengthBuf2 - 1;
     1699      int tmp= tmin (repLengthBuf1, deggSubLg + 1);
     1700      for (int ind= 0; ind < tmp; ind++)
     1701      //{
     1702        nmod_poly_set_coeff_ui (g, ind + lg, n_submod (nmod_poly_get_coeff_ui (g, ind + lg), nmod_poly_get_coeff_ui (buf1, ind), getCharacteristic()));
     1703        //cout << "nmod_poly_get_coeff_ui (g, ind + lg)= " << nmod_poly_get_coeff_ui (g, ind + lg) << "\n";
     1704      //}
     1705    }
     1706    if (lg < 0)
     1707    {
     1708      nmod_poly_clear (buf1);
     1709      nmod_poly_clear (buf2);
     1710      nmod_poly_clear (buf3);
     1711      break;
     1712    }
     1713    if (degfSubLf >= 0)
     1714    {
     1715      for (int ind= 0; ind < repLengthBuf2; ind++)
     1716        nmod_poly_set_coeff_ui (f, ind + lf, n_submod (nmod_poly_get_coeff_ui (f, ind + lf), nmod_poly_get_coeff_ui (buf2, ind), getCharacteristic()));
     1717    }
     1718    nmod_poly_clear (buf1);
     1719    nmod_poly_clear (buf2);
     1720    nmod_poly_clear (buf3);
     1721    /*nmod_poly_init_preinv (buf1, getCharacteristic(), ninv);
     1722    nmod_poly_init_preinv (buf2, getCharacteristic(), ninv);
     1723    nmod_poly_init_preinv (buf3, getCharacteristic(), ninv);*/
     1724  }
     1725
     1726  /*nmod_poly_clear (buf1);
     1727  nmod_poly_clear (buf2);
     1728  nmod_poly_clear (buf3);*/
     1729  nmod_poly_clear (f);
     1730  nmod_poly_clear (g);
     1731
     1732  return result;
     1733}
     1734
     1735CanonicalForm
     1736reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
     1737{
     1738  Variable y= Variable (2);
     1739  Variable x= Variable (1);
     1740
     1741  fmpz_poly_t f, g;
     1742  fmpz_poly_init (f);
     1743  fmpz_poly_init (g);
     1744  fmpz_poly_set (f, F);
     1745  fmpz_poly_set (g, G);
     1746  int degf= fmpz_poly_degree(f);
     1747  int degg= fmpz_poly_degree(g);
     1748
     1749
     1750  fmpz_poly_t buf1,buf2, buf3;
     1751  /*nmod_poly_init_preinv (buf1, getCharacteristic(), ninv);
     1752  nmod_poly_init_preinv (buf2, getCharacteristic(), ninv);
     1753  nmod_poly_init_preinv (buf3, getCharacteristic(), ninv);*/
     1754
     1755  if (fmpz_poly_length (f) < (long) d*(k+1)) //zero padding
     1756    fmpz_poly_fit_length (f,(long)d*(k+1));
     1757
     1758  CanonicalForm result= 0;
     1759  int i= 0;
     1760  int lf= 0;
     1761  int lg= d*k;
     1762  int degfSubLf= degf;
     1763  int deggSubLg= degg-lg;
     1764  int repLengthBuf2;
     1765  int repLengthBuf1;
     1766  //zz_p zzpZero= zz_p();
     1767  fmpz_t tmp1, tmp2;
     1768  while (degf >= lf || lg >= 0)
     1769  {
     1770    if (degfSubLf >= d)
     1771      repLengthBuf1= d;
     1772    else if (degfSubLf < 0)
     1773      repLengthBuf1= 0;
     1774    else
     1775      repLengthBuf1= degfSubLf + 1;
     1776    //cout << "repLengthBuf1= " << repLengthBuf1 << "\n";
     1777    //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
     1778    fmpz_poly_init2 (buf1, repLengthBuf1);
     1779    //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
     1780
     1781    for (int ind= 0; ind < repLengthBuf1; ind++)
     1782    {
     1783      fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
     1784      fmpz_poly_set_coeff_fmpz (buf1, ind, tmp1);
     1785    }
     1786    _fmpz_poly_normalise (buf1);
     1787    //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
     1788
     1789    repLengthBuf1= fmpz_poly_length (buf1);
     1790
     1791    if (deggSubLg >= d - 1)
     1792      repLengthBuf2= d - 1;
     1793    else if (deggSubLg < 0)
     1794      repLengthBuf2= 0;
     1795    else
     1796      repLengthBuf2= deggSubLg + 1;
     1797
     1798    //cout << "repLengthBuf2= " << repLengthBuf2 << "\n";
     1799    //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
     1800    fmpz_poly_init2 (buf2, repLengthBuf2);
     1801    //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
     1802    for (int ind= 0; ind < repLengthBuf2; ind++)
     1803    {
     1804      fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
     1805      fmpz_poly_set_coeff_fmpz (buf2, ind, tmp1);
     1806    }
     1807
     1808    _fmpz_poly_normalise (buf2);
     1809    //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
     1810    repLengthBuf2= fmpz_poly_length (buf2);
     1811
     1812    fmpz_poly_init2 (buf3, repLengthBuf2 + d);
     1813    for (int ind= 0; ind < repLengthBuf1; ind++)
     1814    {
     1815      fmpz_poly_get_coeff_fmpz (tmp1, buf1, ind);  //oder fmpz_set (fmpz_poly_get_coeff_ptr (buf3, ind),fmpz_poly_get_coeff_ptr (buf1, ind))
     1816      fmpz_poly_set_coeff_fmpz (buf3, ind, tmp1);
     1817    }
     1818    for (int ind= repLengthBuf1; ind < d; ind++)
     1819      fmpz_poly_set_coeff_ui (buf3, ind, 0);
     1820    for (int ind= 0; ind < repLengthBuf2; ind++)
     1821    {
     1822      fmpz_poly_get_coeff_fmpz (tmp1, buf2, ind);
     1823      fmpz_poly_set_coeff_fmpz (buf3, ind + d, tmp1);
     1824    }
     1825    _fmpz_poly_normalise (buf3);
     1826
     1827    result += convertFmpz_poly_t2FacCF (buf3, x)*power (y, i);
     1828    i++;
     1829
     1830
     1831    lf= i*d;
     1832    degfSubLf= degf - lf;
     1833
     1834    lg= d*(k-i);
     1835    deggSubLg= degg - lg;
     1836
     1837    if (lg >= 0 && deggSubLg > 0)
     1838    {
     1839      if (repLengthBuf2 > degfSubLf + 1)
     1840        degfSubLf= repLengthBuf2 - 1;
     1841      int tmp= tmin (repLengthBuf1, deggSubLg + 1);
     1842      for (int ind= 0; ind < tmp; ind++)
     1843      {
     1844        fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
     1845        fmpz_poly_get_coeff_fmpz (tmp2, buf1, ind);
     1846        fmpz_sub (tmp1, tmp1, tmp2);
     1847        fmpz_poly_set_coeff_fmpz (g, ind + lg, tmp1);
     1848        //cout << "nmod_poly_get_coeff_ui (g, ind + lg)= " << nmod_poly_get_coeff_ui (g, ind + lg) << "\n";
     1849      }
     1850    }
     1851    if (lg < 0)
     1852    {
     1853      fmpz_poly_clear (buf1);
     1854      fmpz_poly_clear (buf2);
     1855      fmpz_poly_clear (buf3);
     1856      break;
     1857    }
     1858    if (degfSubLf >= 0)
     1859    {
     1860      for (int ind= 0; ind < repLengthBuf2; ind++)
     1861      {
     1862        fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
     1863        fmpz_poly_get_coeff_fmpz (tmp2, buf2, ind);
     1864        fmpz_sub (tmp1, tmp1, tmp2);
     1865        fmpz_poly_set_coeff_fmpz (f, ind + lf, tmp1);
     1866      }
     1867    }
     1868    fmpz_poly_clear (buf1);
     1869    fmpz_poly_clear (buf2);
     1870    fmpz_poly_clear (buf3);
     1871    /*nmod_poly_init_preinv (buf1, getCharacteristic(), ninv);
     1872    nmod_poly_init_preinv (buf2, getCharacteristic(), ninv);
     1873    nmod_poly_init_preinv (buf3, getCharacteristic(), ninv);*/
     1874  }
     1875
     1876  /*nmod_poly_clear (buf1);
     1877  nmod_poly_clear (buf2);
     1878  nmod_poly_clear (buf3);*/
     1879  fmpz_poly_clear (f);
     1880  fmpz_poly_clear (g);
     1881  fmpz_clear (tmp1);
     1882  fmpz_clear (tmp2);
     1883
     1884  return result;
     1885}
     1886
     1887CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
     1888{
     1889  Variable y= Variable (2);
     1890  Variable x= Variable (1);
     1891
     1892  nmod_poly_t f;
     1893  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
     1894  nmod_poly_init_preinv (f, getCharacteristic(), ninv);
     1895  nmod_poly_set (f, F);
     1896
     1897  nmod_poly_t buf;
     1898  //nmod_poly_init_preinv (buf, getCharacteristic(), ninv);
     1899  CanonicalForm result= 0;
     1900  int i= 0;
     1901  int degf= nmod_poly_degree(f);
     1902  int k= 0;
     1903  int degfSubK;
     1904  int repLength;
     1905  while (degf >= k)
     1906  {
     1907    degfSubK= degf - k;
     1908    if (degfSubK >= d)
     1909      repLength= d;
     1910    else
     1911      repLength= degfSubK + 1;
     1912
     1913    nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
     1914    for (int j= 0; j < repLength; j++)
     1915      nmod_poly_set_coeff_ui (buf, j, nmod_poly_get_coeff_ui (f, j + k));
     1916    //printf ("after set\n");
     1917    _nmod_poly_normalise (buf);
     1918    //printf ("after normalise\n");
     1919
     1920    result += convertnmod_poly_t2FacCF (buf, x)*power (y, i);
     1921    i++;
     1922    k= d*i;
     1923    nmod_poly_clear (buf);
     1924    //nmod_poly_init_preinv (buf, getCharacteristic(), ninv);
     1925  }
     1926  //printf ("after while\n");
     1927  //nmod_poly_clear (buf);
     1928  nmod_poly_clear (f);
     1929  //printf ("after clear\n");
     1930
     1931  return result;
     1932}
     1933
     1934CanonicalForm
     1935mulMod2FLINTFpReci (const CanonicalForm& F, const CanonicalForm& G, const
     1936                    CanonicalForm& M)
     1937{
     1938  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
     1939  d1 /= 2;
     1940  d1 += 1;
     1941
     1942  nmod_poly_t F1, F2;
     1943  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
     1944  nmod_poly_init_preinv (F1, getCharacteristic(), ninv);
     1945  nmod_poly_init_preinv (F2, getCharacteristic(), ninv);
     1946  kronSubRecipro (F1, F2, F, d1);
     1947  /*zz_pX TF1, TF2;
     1948  kronSubRecipro (TF1, TF2, F, d1);
     1949  Variable x= Variable (1);
     1950  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF1, x) == convertnmod_poly_t2FacCF (F1, x)) << "\n";
     1951  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF2, x) == convertnmod_poly_t2FacCF (F2, x)) << "\n";*/
     1952
     1953  nmod_poly_t G1, G2;
     1954  nmod_poly_init_preinv (G1, getCharacteristic(), ninv);
     1955  nmod_poly_init_preinv (G2, getCharacteristic(), ninv);
     1956  kronSubRecipro (G1, G2, G, d1);
     1957  /*zz_pX TG1, TG2;
     1958  kronSubRecipro (TG1, TG2, G, d1);
     1959  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TG1, x) == convertnmod_poly_t2FacCF (G1, x)) << "\n";
     1960  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TG2, x) == convertnmod_poly_t2FacCF (G2, x)) << "\n";*/
     1961
     1962
     1963  int k= d1*degree (M);
     1964  nmod_poly_mullow (F1, F1, G1, (long) k);
     1965  /*MulTrunc (TF1, TF1, TG1, (long) k);
     1966  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF1, x) == convertnmod_poly_t2FacCF (F1, x)) << "\n";*/
     1967
     1968  int degtailF= degree (tailcoeff (F), 1);;
     1969  int degtailG= degree (tailcoeff (G), 1);
     1970  int taildegF= taildegree (F);
     1971  int taildegG= taildegree (G);
     1972
     1973  int b= nmod_poly_degree (F2) + nmod_poly_degree (G2) - k - degtailF - degtailG + d1*(2+taildegF + taildegG);
     1974  nmod_poly_mulhigh (F2, F2, G2, b);
     1975  nmod_poly_shift_right (F2, F2, b);
     1976  /*mul (TF2, TF2, TG2);
     1977  if (deg (TF2) > k)
     1978  {
     1979    int b= deg (TF2) - k + 2;
     1980    TF2 >>= b;
     1981  }
     1982  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF2, x) == convertnmod_poly_t2FacCF (F2, x)) << "\n";*/
     1983  int d2= tmax (nmod_poly_degree (F2)/d1, nmod_poly_degree (F1)/d1);
     1984
     1985
     1986  CanonicalForm result= reverseSubst (F1, F2, d1, d2);
     1987  //CanonicalForm NTLres= reverseSubst (TF1, TF2, d1, d2);
     1988  //cout << "FINALtestkronSubReci= " << (NTLres == result) << "\n";
     1989
     1990  nmod_poly_clear (F1);
     1991  nmod_poly_clear (F2);
     1992  nmod_poly_clear (G1);
     1993  nmod_poly_clear (G2);
     1994  return result;
     1995}
     1996
     1997CanonicalForm
     1998mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
     1999                CanonicalForm& M)
     2000{
     2001  CanonicalForm A= F;
     2002  CanonicalForm B= G;
     2003
     2004  int degAx= degree (A, 1);
     2005  int degAy= degree (A, 2);
     2006  int degBx= degree (B, 1);
     2007  int degBy= degree (B, 2);
     2008  int d1= degAx + 1 + degBx;
     2009  int d2= tmax (degAy, degBy);
     2010
     2011  if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
     2012    return mulMod2FLINTFpReci (A, B, M);
     2013
     2014  nmod_poly_t FLINTA, FLINTB;
     2015  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
     2016  nmod_poly_init_preinv (FLINTA, getCharacteristic(), ninv);
     2017  nmod_poly_init_preinv (FLINTB, getCharacteristic(), ninv);
     2018  kronSubFp (FLINTA, A, d1);
     2019  kronSubFp (FLINTB, B, d1);
     2020
     2021  int k= d1*degree (M);
     2022  nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
     2023
     2024  A= reverseSubstFp (FLINTA, d1);
     2025
     2026  nmod_poly_clear (FLINTA);
     2027  nmod_poly_clear (FLINTB);
     2028  return A;
     2029}
     2030
     2031CanonicalForm
     2032mulMod2FLINTQReci (const CanonicalForm& F, const CanonicalForm& G, const
     2033                    CanonicalForm& M)
     2034{
     2035  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
     2036  d1 /= 2;
     2037  d1 += 1;
     2038
     2039  fmpz_poly_t F1, F2;
     2040  fmpz_poly_init (F1);
     2041  fmpz_poly_init (F2);
     2042  kronSubRecipro (F1, F2, F, d1);
     2043  /*zz_pX TF1, TF2;
     2044  kronSubRecipro (TF1, TF2, F, d1);
     2045  Variable x= Variable (1);
     2046  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF1, x) == convertnmod_poly_t2FacCF (F1, x)) << "\n";
     2047  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF2, x) == convertnmod_poly_t2FacCF (F2, x)) << "\n";*/
     2048
     2049  fmpz_poly_t G1, G2;
     2050  fmpz_poly_init (G1);
     2051  fmpz_poly_init (G2);
     2052  kronSubRecipro (G1, G2, G, d1);
     2053  /*zz_pX TG1, TG2;
     2054  kronSubRecipro (TG1, TG2, G, d1);
     2055  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TG1, x) == convertnmod_poly_t2FacCF (G1, x)) << "\n";
     2056  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TG2, x) == convertnmod_poly_t2FacCF (G2, x)) << "\n";*/
     2057
     2058
     2059  int k= d1*degree (M);
     2060  fmpz_poly_mullow (F1, F1, G1, (long) k);
     2061  /*MulTrunc (TF1, TF1, TG1, (long) k);
     2062  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF1, x) == convertnmod_poly_t2FacCF (F1, x)) << "\n";*/
     2063
     2064  int degtailF= degree (tailcoeff (F), 1);;
     2065  int degtailG= degree (tailcoeff (G), 1);
     2066  int taildegF= taildegree (F);
     2067  int taildegG= taildegree (G);
     2068
     2069  int b= fmpz_poly_degree (F2) + fmpz_poly_degree (G2) - k - degtailF - degtailG + d1*(2+taildegF + taildegG);
     2070  fmpz_poly_mulhigh_n (F2, F2, G2, b);
     2071  fmpz_poly_shift_right (F2, F2, b);
     2072  /*mul (TF2, TF2, TG2);
     2073  if (deg (TF2) > k)
     2074  {
     2075    int b= deg (TF2) - k + 2;
     2076    TF2 >>= b;
     2077  }
     2078  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF2, x) == convertnmod_poly_t2FacCF (F2, x)) << "\n";*/
     2079  int d2= tmax (fmpz_poly_degree (F2)/d1, fmpz_poly_degree (F1)/d1);
     2080
     2081
     2082  CanonicalForm result= reverseSubst (F1, F2, d1, d2);
     2083
     2084  //CanonicalForm NTLres= reverseSubst (TF1, TF2, d1, d2);
     2085  //cout << "FINALtestkronSubReci= " << (NTLres == result) << "\n";
     2086
     2087  fmpz_poly_clear (F1);
     2088  fmpz_poly_clear (F2);
     2089  fmpz_poly_clear (G1);
     2090  fmpz_poly_clear (G2);
     2091  return result;
     2092}
     2093
     2094CanonicalForm
     2095mulMod2FLINTQ (const CanonicalForm& F, const CanonicalForm& G, const
     2096                CanonicalForm& M)
     2097{
     2098  CanonicalForm A= F;
     2099  CanonicalForm B= G;
     2100
     2101  int degAx= degree (A, 1);
     2102  //int degAy= degree (A, 2);
     2103  int degBx= degree (B, 1);
     2104  //int degBy= degree (B, 2);
     2105  int d1= degAx + 1 + degBx;
     2106  //int d2= tmax (degAy, degBy);
     2107
     2108  CanonicalForm f= bCommonDen (F);
     2109  CanonicalForm g= bCommonDen (G);
     2110  A *= f;
     2111  B *= g;
     2112
     2113  fmpz_poly_t FLINTA, FLINTB;
     2114  fmpz_poly_init (FLINTA);
     2115  fmpz_poly_init (FLINTB);
     2116  kronSub (FLINTA, A, d1);
     2117  kronSub (FLINTB, B, d1);
     2118  int k= d1*degree (M);
     2119
     2120  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
     2121  A= reverseSubstQ (FLINTA, d1);
     2122  fmpz_poly_clear (FLINTA);
     2123  fmpz_poly_clear (FLINTB);
     2124  return A/(f*g);
     2125}
     2126
     2127#endif
    14532128
    14542129CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
     
    14862161  if (sizeF < fallBackToNaive || sizeG < fallBackToNaive)
    14872162    return mod (F*G, M);
     2163
     2164#ifdef HAVE_FLINT
     2165  Variable alpha;
     2166  if (getCharacteristic() == 0 && !hasFirstAlgVar (F, alpha) && ! hasFirstAlgVar (G, alpha))
     2167  {
     2168    CanonicalForm result= mulMod2FLINTQ (F, G, M);
     2169    CanonicalForm test= mod (F*G, M);
     2170    if (test != result)
     2171      printf ("GOTTVERDAMMT\n");
     2172    return result;
     2173    //return mulMod2FLINTQ (F, G, M);
     2174  }
     2175#endif
    14882176
    14892177  if (getCharacteristic() > 0 && CFFactory::gettype() != GaloisFieldDomain &&
Note: See TracChangeset for help on using the changeset viewer.