Changeset 81d96c in git


Ignore:
Timestamp:
Feb 9, 2012, 5:59:15 PM (11 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
Children:
c0851051c3f2502093559a89e463f364c0816aa0
Parents:
0e2e23228417d64209e60f3d18fa1bd11233839d
git-author:
Martin Lee <martinlee84@web.de>2012-02-09 17:59:15+01:00
git-committer:
Martin Lee <martinlee84@web.de>2012-04-04 14:42:25+02:00
Message:
chg: more reorganization
Location:
factory
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • factory/cf_gcd_smallp.cc

    r0e2e23 r81d96c  
    2626
    2727#include "canonicalform.h"
     28#include "algext.h"
    2829#include "cf_map.h"
    2930#include "cf_util.h"
    … …  
    40324033  lcs [0]= shiftedLCsEval1.getFirst();
    40334034  lcs [1]= shiftedLCsEval2.getFirst();
    4034   henselLift122 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
    4035                  lcs, false);
     4035  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
     4036                        lcs, false);
    40364037
    40374038  for (CFListIterator i= factors; i.hasItem(); i++)
    … …  
    40494050    for (int i= 1; i < U.level() - 1; i++)
    40504051      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
    4051     factors= henselLift2 (UEval, factors, liftBounds, U.level() - 1, false,
    4052                           shiftedLCsEval1, shiftedLCsEval2, Pi, diophant,
    4053                           noOneToOne);
     4052    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
     4053                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
     4054                                  diophant, noOneToOne);
    40544055    delete [] liftBounds;
    40554056    if (noOneToOne)

  • factory/facFactorize.cc

    r0e2e23 r81d96c  
    481481      CFArray Pi;
    482482      CFList diophant;
    483       henselLift122 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
    484                      leadingCoeffs, false);
     483      nonMonicHenselLift12 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
     484                            leadingCoeffs, false);
    485485
    486486      if (sqrfPartF.level() > 2)

  • factory/facFqBivarUtil.cc

    r0e2e23 r81d96c  
    1616
    1717#include "cf_map.h"
     18#include "algext.h"
    1819#include "cf_map_ext.h"
    1920#include "templates/ftmpl_functions.h"

  • factory/facFqFactorize.cc

    r0e2e23 r81d96c  
    16581658      CFArray Pi;
    16591659      CFList diophant;
    1660       henselLift122 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
    1661                      leadingCoeffs, false);
     1660      nonMonicHenselLift12 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
     1661                            leadingCoeffs, false);
    16621662
    16631663      if (sqrfPartF.level() > 2)

  • factory/facHensel.cc

    r0e2e23 r81d96c  
    2525#include "timing.h"
    2626
     27#include "algext.h"
    2728#include "facHensel.h"
    2829#include "facMul.h"
    29 #include "cf_util.h"
    3030#include "fac_util.h"
    3131#include "cf_algorithm.h"
    … …  
    125125}
    126126
    127 static inline
     127static
    128128CFList mapinto (const CFList& L)
    129129{
    … …  
    134134}
    135135
    136 static inline
     136static
    137137int mod (const CFList& L, const CanonicalForm& p)
    138138{
    … …  
    146146
    147147
    148 static inline void
     148static void
    149149chineseRemainder (const CFList & x1, const CanonicalForm & q1,
    150150                  const CFList & x2, const CanonicalForm & q2,
    … …  
    162162}
    163163
    164 static inline
     164static
    165165CFList Farey (const CFList& L, const CanonicalForm& q)
    166166{
    … …  
    171171}
    172172
    173 static inline
     173static
    174174CFList replacevar (const CFList& L, const Variable& a, const Variable& b)
    175175{
    … …  
    359359}
    360360
    361 static inline
    362361CFList diophantine (const CanonicalForm& F, const CFList& factors)
    363362{
    … …  
    669668}
    670669
    671 static inline
    672670CFList
    673 biDiophantine (const CanonicalForm& F, const CFList& factors, const int d)
     671biDiophantine (const CanonicalForm& F, const CFList& factors, int d)
    674672{
    675673  Variable y= F.mvar();
    … …  
    770768}
    771769
    772 static inline
    773770CFList
    774771multiRecDiophantine (const CanonicalForm& F, const CFList& factors,
    775                      const CFList& recResult, const CFList& M, const int d)
     772                     const CFList& recResult, const CFList& M, int d)
    776773{
    777774  Variable y= F.mvar();
    … …  
    859856}
    860857
    861 static inline
     858static
    862859void
    863860henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
    … …  
    11011098
    11021099CFList
    1103 henselLift23 (const CFList& eval, const CFList& factors, const int* l, CFList&
     1100henselLift23 (const CFList& eval, const CFList& factors, int* l, CFList&
    11041101              diophant, CFArray& Pi, CFMatrix& M)
    11051102{
    … …  
    11681165CFList
    11691166henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
    1170             diophant, CFArray& Pi, CFMatrix& M, const int lOld, const int
    1171             lNew)
     1167            diophant, CFArray& Pi, CFMatrix& M, int lOld, int lNew)
    11721168{
    11731169  diophant= multiRecDiophantine (F.getFirst(), factors, diophant, MOD, lOld);
    … …  
    12091205
    12101206CFList
    1211 henselLift (const CFList& eval, const CFList& factors, const int* l, const int
    1212             lLength, bool sort)
     1207henselLift (const CFList& eval, const CFList& factors, int* l, int lLength,
     1208            bool sort)
    12131209{
    12141210  CFList diophant;
    … …  
    12431239}
    12441240
     1241// nonmonic
     1242
    12451243void
    1246 henselStep122 (const CanonicalForm& F, const CFList& factors,
    1247               CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
    1248               CFArray& Pi, int j, const CFArray& /*LCs*/)
     1244nonMonicHenselStep12 (const CanonicalForm& F, const CFList& factors,
     1245                      CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
     1246                      CFArray& Pi, int j, const CFArray& /*LCs*/)
    12491247{
    12501248  Variable x= F.mvar();
    … …  
    14551453
    14561454void
    1457 henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
    1458               CFList& diophant, CFMatrix& M, const CFArray& LCs, bool sort)
     1455nonMonicHenselLift12 (const CanonicalForm& F, CFList& factors, int l,
     1456                      CFArray& Pi, CFList& diophant, CFMatrix& M,
     1457                      const CFArray& LCs, bool sort)
    14591458{
    14601459  if (sort)
    … …  
    15221521
    15231522  for (i= 1; i < l; i++)
    1524     henselStep122 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
     1523    nonMonicHenselStep12 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
    15251524
    15261525  factors= CFList();
    … …  
    15331532/// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$
    15341533/// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$
    1535 static inline
    15361534CFList
    15371535diophantine (const CFList& recResult, const CFList& factors,
    … …  
    16121610}
    16131611
    1614 static inline
    16151612void
    1616 henselStep2 (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
    1617             const CFList& diophant, CFMatrix& M, CFArray& Pi,
    1618             const CFList& products, int j, const CFList& MOD, bool& noOneToOne)
     1613nonMonicHenselStep (const CanonicalForm& F, const CFList& factors,
     1614                    CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
     1615                    CFArray& Pi, const CFList& products, int j,
     1616                    const CFList& MOD, bool& noOneToOne)
    16191617{
    16201618  CanonicalForm E;
    … …  
    18521850
    18531851CFList
    1854 henselLift232 (const CFList& eval, const CFList& factors, int* l, CFList&
    1855               diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
    1856               const CFList& LCs2, bool& bad)
     1852nonMonicHenselLift232(const CFList& eval, const CFList& factors, int* l, CFList&
     1853                      diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
     1854                      const CFList& LCs2, bool& bad)
    18571855{
    18581856  CFList buf= factors;
    … …  
    19021900  for (int d= 1; d < l[1]; d++)
    19031901  {
    1904     henselStep2 (j.getItem(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
     1902    nonMonicHenselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, products,
     1903                        d, MOD, bad);
    19051904    if (bad)
    19061905      return CFList();
    … …  
    19141913
    19151914CFList
    1916 henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
    1917             diophant, CFArray& Pi, CFMatrix& M, const int lOld, int&
    1918             lNew, const CFList& LCs1, const CFList& LCs2, bool& bad)
     1915nonMonicHenselLift2 (const CFList& F, const CFList& factors, const CFList& MOD,
     1916                    CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
     1917                    int& lNew, const CFList& LCs1, const CFList& LCs2, bool& bad
     1918                   )
    19191919{
    19201920  int k= 0;
    … …  
    19631963  for (int d= 1; d < lNew; d++)
    19641964  {
    1965     henselStep2 (F.getLast(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
     1965    nonMonicHenselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, products,
     1966                        d, MOD, bad);
    19661967    if (bad)
    19671968      return CFList();
    … …  
    19751976
    19761977CFList
    1977 henselLift2 (const CFList& eval, const CFList& factors, int* l, const int
    1978              lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
    1979              const CFArray& Pi, const CFList& diophant, bool& bad)
     1978nonMonicHenselLift2 (const CFList& eval, const CFList& factors, int* l, int
     1979                    lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
     1980                    const CFArray& Pi, const CFList& diophant, bool& bad)
    19801981{
    19811982  CFList bufDiophant= diophant;
    … …  
    19851986  CFArray bufPi= Pi;
    19861987  CFMatrix M= CFMatrix (l[1], factors.length());
    1987   CFList result= henselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2,
    1988                               bad);
     1988  CFList result=
     1989    nonMonicHenselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2, bad);
    19891990  if (bad)
    19901991    return CFList();
    … …  
    20142015    bufLCs2.append (jjj.getItem());
    20152016    M= CFMatrix (l[i], factors.length());
    2016     result= henselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M, l[i - 1],
    2017                          l[i], bufLCs1, bufLCs2, bad);
     2017    result= nonMonicHenselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M,
     2018                                 l[i - 1], l[i], bufLCs1, bufLCs2, bad);
    20182019    if (bad)
    20192020      return CFList();
    … …  
    21122113  for (int d= 1; d < liftBound; d++)
    21132114  {
    2114     henselStep2 (F, factors, bufFactors, diophant, M, Pi, products, d, MOD, bad);
     2115    nonMonicHenselStep (F, factors, bufFactors, diophant, M, Pi, products, d,
     2116                        MOD, bad);
    21152117    if (bad)
    21162118      return CFList();
    … …  
    21252127CFList
    21262128nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs,
    2127                     CFList& diophant, CFArray& Pi, CFMatrix& M, const int lOld,
     2129                    CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
    21282130                    int& lNew, const CFList& MOD, bool& noOneToOne
    21292131                   )
    … …  
    21922194  for (int d= 1; d < lNew; d++)
    21932195  {
    2194     henselStep2 (F.getLast(), factors, bufFactors, diophant, M, Pi, products, d,
    2195                  MOD, noOneToOne);
     2196    nonMonicHenselStep (F.getLast(), factors, bufFactors, diophant, M, Pi,
     2197                        products, d, MOD, noOneToOne);
    21962198    if (noOneToOne)
    21972199      return CFList();

  • factory/facHensel.h

    r0e2e23 r81d96c  
    2525
    2626#include "canonicalform.h"
    27 #include "cf_iter.h"
    28 #include "templates/ftmpl_functions.h"
    29 #include "algext.h"
    3027
    3128#ifdef HAVE_NTL
    … …  
    8784                                     ///< including leading coefficient
    8885                                     ///< as first element.
    89               const int* l,          ///< [in] l[0]: precision of bivariate
     86              int* l,          ///< [in] l[0]: precision of bivariate
    9087                                     ///< lifting, l[1] as above
    9188              CFList& diophant,      ///< [in,out] returns the result of
    … …  
    131128            CFArray& Pi,           ///< [in,out] stores intermediate results
    132129            CFMatrix& M,           ///< [in,out] stores intermediate results
    133             const int lOld,       ///< [in] lifting precision of F.mvar()
    134             const int lNew        ///< [in] lifting precision of G.mvar()
     130            int lOld,       ///< [in] lifting precision of F.mvar()
     131            int lNew        ///< [in] lifting precision of G.mvar()
    135132           );
    136133
    … …  
    148145            const CFList& factors, ///< [in] bivariate factors, including
    149146                                   ///< leading coefficient
    150             const int* l,          ///< [in] lifting bounds
    151             const int lLength,     ///< [in] length of l
     147            int* l,          ///< [in] lifting bounds
     148            int lLength,     ///< [in] length of l
    152149            bool sort= true        ///< [in] sort factors by degree in
    153150                                   ///< Variable(1)
    154151           );
    155152
    156 /// two factor Hensel lifting from univariate to bivariate, factors need not to
    157 /// be monic
    158 void
    159 henselLift122 (const CanonicalForm& F,///< [in] a bivariate poly
     153/// Hensel lifting from univariate to bivariate, factors need not to be monic
     154void
     155nonMonicHenselLift12 (const CanonicalForm& F,///< [in] a bivariate poly
    160156               CFList& factors,       ///< [in, out] a list of univariate polys
    161157                                      ///< lifted factors
    … …  
    174170/// @return @a henselLift122 returns a list of lifted factors
    175171CFList
    176 henselLift2 (const CFList& eval,   ///< [in] a list of polynomials the last
     172nonMonicHenselLift2 (const CFList& eval,///< [in] a list of polynomials the last
    177173                                   ///< element is a compressed multivariate
    178174                                   ///< poly, last but one element equals the
    … …  
    182178             const CFList& factors,///< [in] bivariate factors
    183179             int* l,               ///< [in] lift bounds
    184              const int lLength,    ///< [in] length of l
     180             int lLength,          ///< [in] length of l
    185181             bool sort,            ///< [in] if true factors are sorted by
    186182                                   ///< their degree in Variable(1)

  • factory/facMul.cc

    r0e2e23 r81d96c  
    1717#include "facMul.h"
    1818#include "algext.h"
     19#include "cf_util.h"
    1920#include "templates/ftmpl_functions.h"
    2021
    … …  
    2627#include "FLINTconvert.h"
    2728#endif
     29
     30// univariate polys
    2831
    2932#ifdef HAVE_FLINT
    … …  
    4851
    4952CanonicalForm
    50 reverseSubst (const fmpz_poly_t F, int d, const Variable& alpha,
    51               const CanonicalForm& den)
     53reverseSubstQa (const fmpz_poly_t F, int d, const Variable& alpha,
     54                const CanonicalForm& den)
    5255{
    5356  Variable x= Variable (1);
    … …  
    111114
    112115  denA *= denB;
    113   A= reverseSubst (FLINTA, d, alpha, denA);
     116  A= reverseSubstQa (FLINTA, d, alpha, denA);
    114117
    115118  fmpz_poly_clear (FLINTA);
    … …  
    222225
    223226  denA *= denB;
    224   A= reverseSubst (FLINTA, d, alpha, denA);
     227  A= reverseSubstQa (FLINTA, d, alpha, denA);
    225228  fmpz_poly_clear (FLINTA);
    226229  fmpz_poly_clear (FLINTB);
    … …  
    570573}
    571574
     575// end univariate polys
     576//*************************
     577// bivariate polys
     578
    572579#ifdef HAVE_FLINT
    573580void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
    … …  
    592599}
    593600
    594 /*void kronSubQ (fmpz_poly_t result, const CanonicalForm& A, int d)
    595 {
    596   int degAy= degree (A);
    597   fmpz_poly_init2 (result, d*(degAy + 1));
    598   _fmpz_poly_set_length (result, d*(degAy+1));
    599 
    600   CFIterator j;
    601   for (CFIterator i= A; i.hasTerms(); i++)
    602   {
    603     if (i.coeff().inBas
    604     convertFacCF2Fmpz_poly_t (buf, i.coeff());
    605 
    606     int k= i.exp()*d;
    607     int bufRepLength= (int) fmpz_poly_length (buf);
    608     for (int j= 0; j < bufRepLength; j++)
    609     {
    610       fmpz_poly_get_coeff_fmpz (coeff, buf, j);
    611       fmpz_poly_set_coeff_fmpz (result, j + k, coeff);
    612     }
    613     fmpz_poly_clear (buf);
    614   }
    615   fmpz_clear (coeff);
    616   _fmpz_poly_normalise (result);
    617 }*/
    618 
    619 // A is a bivariate poly over Qa!!!!
    620 // d2= 2*deg_alpha + 1
    621 // d1= 2*deg_x*d2+1
    622601void kronSubQa (fmpq_poly_t result, const CanonicalForm& A, int d1, int d2)
    623602{
    … …  
    664643  _fmpq_poly_normalise (result);
    665644}
    666 #endif
    667 
    668 zz_pX kronSubFp (const CanonicalForm& A, int d)
    669 {
    670   int degAy= degree (A);
    671   zz_pX result;
    672   result.rep.SetLength (d*(degAy + 1));
    673 
    674   zz_p *resultp;
    675   resultp= result.rep.elts();
    676   zz_pX buf;
    677   zz_p *bufp;
    678   int j, k, bufRepLength;
    679 
    680   for (CFIterator i= A; i.hasTerms(); i++)
    681   {
    682     if (i.coeff().inCoeffDomain())
    683       buf= convertFacCF2NTLzzpX (i.coeff());
    684     else
    685       buf= convertFacCF2NTLzzpX (i.coeff());
    686 
    687     k= i.exp()*d;
    688     bufp= buf.rep.elts();
    689     bufRepLength= (int) buf.rep.length();
    690     for (j= 0; j < bufRepLength; j++)
    691       resultp [j + k]= bufp [j];
    692   }
    693   result.normalize();
    694 
    695   return result;
    696 }
    697 
    698 zz_pEX kronSub (const CanonicalForm& A, int d, const Variable& alpha)
    699 {
    700   int degAy= degree (A);
    701   zz_pEX result;
    702   result.rep.SetLength (d*(degAy + 1));
    703 
    704   Variable v;
    705   zz_pE *resultp;
    706   resultp= result.rep.elts();
    707   zz_pEX buf1;
    708   zz_pE *buf1p;
    709   zz_pX buf2;
    710   zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
    711   int j, k, buf1RepLength;
    712 
    713   for (CFIterator i= A; i.hasTerms(); i++)
    714   {
    715     if (i.coeff().inCoeffDomain())
    716     {
    717       buf2= convertFacCF2NTLzzpX (i.coeff());
    718       buf1= to_zz_pEX (to_zz_pE (buf2));
    719     }
    720     else
    721       buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
    722 
    723     k= i.exp()*d;
    724     buf1p= buf1.rep.elts();
    725     buf1RepLength= (int) buf1.rep.length();
    726     for (j= 0; j < buf1RepLength; j++)
    727       resultp [j + k]= buf1p [j];
    728   }
    729   result.normalize();
    730 
    731   return result;
    732 }
    733645
    734646void
    735 kronSubRecipro (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
    736                 const Variable& alpha)
    737 {
    738   int degAy= degree (A);
    739   subA1.rep.SetLength ((long) d*(degAy + 2));
    740   subA2.rep.SetLength ((long) d*(degAy + 2));
    741 
    742   Variable v;
    743   zz_pE *subA1p;
    744   zz_pE *subA2p;
    745   subA1p= subA1.rep.elts();
    746   subA2p= subA2.rep.elts();
    747   zz_pEX buf;
    748   zz_pE *bufp;
    749   zz_pX buf2;
    750   zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
    751   int j, k, kk, bufRepLength;
    752 
    753   for (CFIterator i= A; i.hasTerms(); i++)
    754   {
    755     if (i.coeff().inCoeffDomain())
    756     {
    757       buf2= convertFacCF2NTLzzpX (i.coeff());
    758       buf= to_zz_pEX (to_zz_pE (buf2));
    759     }
    760     else
    761       buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
    762 
    763     k= i.exp()*d;
    764     kk= (degAy - i.exp())*d;
    765     bufp= buf.rep.elts();
    766     bufRepLength= (int) buf.rep.length();
    767     for (j= 0; j < bufRepLength; j++)
    768     {
    769       subA1p [j + k] += bufp [j];
    770       subA2p [j + kk] += bufp [j];
    771     }
    772   }
    773   subA1.normalize();
    774   subA2.normalize();
    775 }
    776 
    777 void
    778 kronSubRecipro (zz_pX& subA1, zz_pX& subA2,const CanonicalForm& A, int d)
    779 {
    780   int degAy= degree (A);
    781   subA1.rep.SetLength ((long) d*(degAy + 2));
    782   subA2.rep.SetLength ((long) d*(degAy + 2));
    783 
    784   zz_p *subA1p;
    785   zz_p *subA2p;
    786   subA1p= subA1.rep.elts();
    787   subA2p= subA2.rep.elts();
    788   zz_pX buf;
    789   zz_p *bufp;
    790   int j, k, kk, bufRepLength;
    791 
    792   for (CFIterator i= A; i.hasTerms(); i++)
    793   {
    794     buf= convertFacCF2NTLzzpX (i.coeff());
    795 
    796     k= i.exp()*d;
    797     kk= (degAy - i.exp())*d;
    798     bufp= buf.rep.elts();
    799     bufRepLength= (int) buf.rep.length();
    800     for (j= 0; j < bufRepLength; j++)
    801     {
    802       subA1p [j + k] += bufp [j];
    803       subA2p [j + kk] += bufp [j];
    804     }
    805   }
    806   subA1.normalize();
    807   subA2.normalize();
    808 }
    809 
    810 CanonicalForm
    811 reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
    812               const Variable& alpha)
    813 {
    814   Variable y= Variable (2);
    815   Variable x= Variable (1);
    816 
    817   zz_pEX f= F;
    818   zz_pEX g= G;
    819   int degf= deg(f);
    820   int degg= deg(g);
    821 
    822   zz_pEX buf1;
    823   zz_pEX buf2;
    824   zz_pEX buf3;
    825 
    826   zz_pE *buf1p;
    827   zz_pE *buf2p;
    828   zz_pE *buf3p;
    829   if (f.rep.length() < (long) d*(k+1)) //zero padding
    830     f.rep.SetLength ((long)d*(k+1));
    831 
    832   zz_pE *gp= g.rep.elts();
    833   zz_pE *fp= f.rep.elts();
    834   CanonicalForm result= 0;
    835   int i= 0;
    836   int lf= 0;
    837   int lg= d*k;
    838   int degfSubLf= degf;
    839   int deggSubLg= degg-lg;
    840   int repLengthBuf2, repLengthBuf1, ind, tmp;
    841   zz_pE zzpEZero= zz_pE();
    842 
    843   while (degf >= lf || lg >= 0)
    844   {
    845     if (degfSubLf >= d)
    846       repLengthBuf1= d;
    847     else if (degfSubLf < 0)
    848       repLengthBuf1= 0;
    849     else
    850       repLengthBuf1= degfSubLf + 1;
    851     buf1.rep.SetLength((long) repLengthBuf1);
    852 
    853     buf1p= buf1.rep.elts();
    854     for (ind= 0; ind < repLengthBuf1; ind++)
    855       buf1p [ind]= fp [ind + lf];
    856     buf1.normalize();
    857 
    858     repLengthBuf1= buf1.rep.length();
    859 
    860     if (deggSubLg >= d - 1)
    861       repLengthBuf2= d - 1;
    862     else if (deggSubLg < 0)
    863       repLengthBuf2= 0;
    864     else
    865       repLengthBuf2= deggSubLg + 1;
    866 
    867     buf2.rep.SetLength ((long) repLengthBuf2);
    868     buf2p= buf2.rep.elts();
    869     for (ind= 0; ind < repLengthBuf2; ind++)
    870       buf2p [ind]= gp [ind + lg];
    871     buf2.normalize();
    872 
    873     repLengthBuf2= buf2.rep.length();
    874 
    875     buf3.rep.SetLength((long) repLengthBuf2 + d);
    876     buf3p= buf3.rep.elts();
    877     buf2p= buf2.rep.elts();
    878     buf1p= buf1.rep.elts();
    879     for (ind= 0; ind < repLengthBuf1; ind++)
    880       buf3p [ind]= buf1p [ind];
    881     for (ind= repLengthBuf1; ind < d; ind++)
    882       buf3p [ind]= zzpEZero;
    883     for (ind= 0; ind < repLengthBuf2; ind++)
    884       buf3p [ind + d]= buf2p [ind];
    885     buf3.normalize();
    886 
    887     result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
    888     i++;
    889 
    890 
    891     lf= i*d;
    892     degfSubLf= degf - lf;
    893 
    894     lg= d*(k-i);
    895     deggSubLg= degg - lg;
    896 
    897     buf1p= buf1.rep.elts();
    898 
    899     if (lg >= 0 && deggSubLg > 0)
    900     {
    901       if (repLengthBuf2 > degfSubLf + 1)
    902         degfSubLf= repLengthBuf2 - 1;
    903       tmp= tmin (repLengthBuf1, deggSubLg + 1);
    904       for (ind= 0; ind < tmp; ind++)
    905         gp [ind + lg] -= buf1p [ind];
    906     }
    907 
    908     if (lg < 0)
    909       break;
    910 
    911     buf2p= buf2.rep.elts();
    912     if (degfSubLf >= 0)
    913     {
    914       for (ind= 0; ind < repLengthBuf2; ind++)
    915         fp [ind + lf] -= buf2p [ind];
    916     }
    917   }
    918 
    919   return result;
    920 }
    921 
    922 CanonicalForm
    923 reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
    924 {
    925   Variable y= Variable (2);
    926   Variable x= Variable (1);
    927 
    928   zz_pX f= F;
    929   zz_pX g= G;
    930   int degf= deg(f);
    931   int degg= deg(g);
    932 
    933   zz_pX buf1;
    934   zz_pX buf2;
    935   zz_pX buf3;
    936 
    937   zz_p *buf1p;
    938   zz_p *buf2p;
    939   zz_p *buf3p;
    940 
    941   if (f.rep.length() < (long) d*(k+1)) //zero padding
    942     f.rep.SetLength ((long)d*(k+1));
    943 
    944   zz_p *gp= g.rep.elts();
    945   zz_p *fp= f.rep.elts();
    946   CanonicalForm result= 0;
    947   int i= 0;
    948   int lf= 0;
    949   int lg= d*k;
    950   int degfSubLf= degf;
    951   int deggSubLg= degg-lg;
    952   int repLengthBuf2, repLengthBuf1, ind, tmp;
    953   zz_p zzpZero= zz_p();
    954   while (degf >= lf || lg >= 0)
    955   {
    956     if (degfSubLf >= d)
    957       repLengthBuf1= d;
    958     else if (degfSubLf < 0)
    959       repLengthBuf1= 0;
    960     else
    961       repLengthBuf1= degfSubLf + 1;
    962     buf1.rep.SetLength((long) repLengthBuf1);
    963 
    964     buf1p= buf1.rep.elts();
    965     for (ind= 0; ind < repLengthBuf1; ind++)
    966       buf1p [ind]= fp [ind + lf];
    967     buf1.normalize();
    968 
    969     repLengthBuf1= buf1.rep.length();
    970 
    971     if (deggSubLg >= d - 1)
    972       repLengthBuf2= d - 1;
    973     else if (deggSubLg < 0)
    974       repLengthBuf2= 0;
    975     else
    976       repLengthBuf2= deggSubLg + 1;
    977 
    978     buf2.rep.SetLength ((long) repLengthBuf2);
    979     buf2p= buf2.rep.elts();
    980     for (ind= 0; ind < repLengthBuf2; ind++)
    981       buf2p [ind]= gp [ind + lg];
    982 
    983     buf2.normalize();
    984 
    985     repLengthBuf2= buf2.rep.length();
    986 
    987 
    988     buf3.rep.SetLength((long) repLengthBuf2 + d);
    989     buf3p= buf3.rep.elts();
    990     buf2p= buf2.rep.elts();
    991     buf1p= buf1.rep.elts();
    992     for (ind= 0; ind < repLengthBuf1; ind++)
    993       buf3p [ind]= buf1p [ind];
    994     for (ind= repLengthBuf1; ind < d; ind++)
    995       buf3p [ind]= zzpZero;
    996     for (ind= 0; ind < repLengthBuf2; ind++)
    997       buf3p [ind + d]= buf2p [ind];
    998     buf3.normalize();
    999 
    1000     result += convertNTLzzpX2CF (buf3, x)*power (y, i);
    1001     i++;
    1002 
    1003 
    1004     lf= i*d;
    1005     degfSubLf= degf - lf;
    1006 
    1007     lg= d*(k-i);
    1008     deggSubLg= degg - lg;
    1009 
    1010     buf1p= buf1.rep.elts();
    1011 
    1012     if (lg >= 0 && deggSubLg > 0)
    1013     {
    1014       if (repLengthBuf2 > degfSubLf + 1)
    1015         degfSubLf= repLengthBuf2 - 1;
    1016       tmp= tmin (repLengthBuf1, deggSubLg + 1);
    1017       for (ind= 0; ind < tmp; ind++)
    1018         gp [ind + lg] -= buf1p [ind];
    1019     }
    1020     if (lg < 0)
    1021       break;
    1022 
    1023     buf2p= buf2.rep.elts();
    1024     if (degfSubLf >= 0)
    1025     {
    1026       for (ind= 0; ind < repLengthBuf2; ind++)
    1027         fp [ind + lf] -= buf2p [ind];
    1028     }
    1029   }
    1030 
    1031   return result;
    1032 }
    1033 
    1034 CanonicalForm reverseSubst (const zz_pEX& F, int d, const Variable& alpha)
    1035 {
    1036   Variable y= Variable (2);
    1037   Variable x= Variable (1);
    1038 
    1039   zz_pEX f= F;
    1040   zz_pE *fp= f.rep.elts();
    1041 
    1042   zz_pEX buf;
    1043   zz_pE *bufp;
    1044   CanonicalForm result= 0;
    1045   int i= 0;
    1046   int degf= deg(f);
    1047   int k= 0;
    1048   int degfSubK, repLength, j;
    1049   while (degf >= k)
    1050   {
    1051     degfSubK= degf - k;
    1052     if (degfSubK >= d)
    1053       repLength= d;
    1054     else
    1055       repLength= degfSubK + 1;
    1056 
    1057     buf.rep.SetLength ((long) repLength);
    1058     bufp= buf.rep.elts();
    1059     for (j= 0; j < repLength; j++)
    1060       bufp [j]= fp [j + k];
    1061     buf.normalize();
    1062 
    1063     result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
    1064     i++;
    1065     k= d*i;
    1066   }
    1067 
    1068   return result;
    1069 }
    1070 
    1071 CanonicalForm reverseSubstFp (const zz_pX& F, int d)
    1072 {
    1073   Variable y= Variable (2);
    1074   Variable x= Variable (1);
    1075 
    1076   zz_pX f= F;
    1077   zz_p *fp= f.rep.elts();
    1078 
    1079   zz_pX buf;
    1080   zz_p *bufp;
    1081   CanonicalForm result= 0;
    1082   int i= 0;
    1083   int degf= deg(f);
    1084   int k= 0;
    1085   int degfSubK, repLength, j;
    1086   while (degf >= k)
    1087   {
    1088     degfSubK= degf - k;
    1089     if (degfSubK >= d)
    1090       repLength= d;
    1091     else
    1092       repLength= degfSubK + 1;
    1093 
    1094     buf.rep.SetLength ((long) repLength);
    1095     bufp= buf.rep.elts();
    1096     for (j= 0; j < repLength; j++)
    1097       bufp [j]= fp [j + k];
    1098     buf.normalize();
    1099 
    1100     result += convertNTLzzpX2CF (buf, x)*power (y, i);
    1101     i++;
    1102     k= d*i;
    1103   }
    1104 
    1105   return result;
    1106 }
    1107 
    1108 // assumes input to be reduced mod M and to be an element of Fq not Fp
    1109 CanonicalForm
    1110 mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
    1111                   CanonicalForm& M)
    1112 {
    1113   int d1= degree (F, 1) + degree (G, 1) + 1;
    1114   d1 /= 2;
    1115   d1 += 1;
    1116 
    1117   zz_pX F1, F2;
    1118   kronSubRecipro (F1, F2, F, d1);
    1119   zz_pX G1, G2;
    1120   kronSubRecipro (G1, G2, G, d1);
    1121 
    1122   int k= d1*degree (M);
    1123   MulTrunc (F1, F1, G1, (long) k);
    1124 
    1125   int degtailF= degree (tailcoeff (F), 1);
    1126   int degtailG= degree (tailcoeff (G), 1);
    1127   int taildegF= taildegree (F);
    1128   int taildegG= taildegree (G);
    1129   int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
    1130 
    1131   reverse (F2, F2);
    1132   reverse (G2, G2);
    1133   MulTrunc (F2, F2, G2, b + 1);
    1134   reverse (F2, F2, b);
    1135 
    1136   int d2= tmax (deg (F2)/d1, deg (F1)/d1);
    1137   return reverseSubst (F1, F2, d1, d2);
    1138 }
    1139 
    1140 //Kronecker substitution
    1141 CanonicalForm
    1142 mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
    1143               CanonicalForm& M)
    1144 {
    1145   CanonicalForm A= F;
    1146   CanonicalForm B= G;
    1147 
    1148   int degAx= degree (A, 1);
    1149   int degAy= degree (A, 2);
    1150   int degBx= degree (B, 1);
    1151   int degBy= degree (B, 2);
    1152   int d1= degAx + 1 + degBx;
    1153   int d2= tmax (degAy, degBy);
    1154 
    1155   if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
    1156     return mulMod2NTLFpReci (A, B, M);
    1157 
    1158   zz_pX NTLA= kronSubFp (A, d1);
    1159   zz_pX NTLB= kronSubFp (B, d1);
    1160 
    1161   int k= d1*degree (M);
    1162   MulTrunc (NTLA, NTLA, NTLB, (long) k);
    1163 
    1164   A= reverseSubstFp (NTLA, d1);
    1165 
    1166   return A;
    1167 }
    1168 
    1169 // assumes input to be reduced mod M and to be an element of Fq not Fp
    1170 CanonicalForm
    1171 mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
    1172                   CanonicalForm& M, const Variable& alpha)
    1173 {
    1174   int d1= degree (F, 1) + degree (G, 1) + 1;
    1175   d1 /= 2;
    1176   d1 += 1;
    1177 
    1178   zz_pEX F1, F2;
    1179   kronSubRecipro (F1, F2, F, d1, alpha);
    1180   zz_pEX G1, G2;
    1181   kronSubRecipro (G1, G2, G, d1, alpha);
    1182 
    1183   int k= d1*degree (M);
    1184   MulTrunc (F1, F1, G1, (long) k);
    1185 
    1186   int degtailF= degree (tailcoeff (F), 1);
    1187   int degtailG= degree (tailcoeff (G), 1);
    1188   int taildegF= taildegree (F);
    1189   int taildegG= taildegree (G);
    1190   int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
    1191 
    1192   reverse (F2, F2);
    1193   reverse (G2, G2);
    1194   MulTrunc (F2, F2, G2, b + 1);
    1195   reverse (F2, F2, b);
    1196 
    1197   int d2= tmax (deg (F2)/d1, deg (F1)/d1);
    1198   return reverseSubst (F1, F2, d1, d2, alpha);
    1199 }
    1200 
    1201 #ifdef HAVE_FLINT
    1202 CanonicalForm
    1203 mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
    1204                 CanonicalForm& M);
    1205 #endif
    1206 
    1207 CanonicalForm
    1208 mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
    1209               CanonicalForm& M)
    1210 {
    1211   Variable alpha;
    1212   CanonicalForm A= F;
    1213   CanonicalForm B= G;
    1214 
    1215   if (hasFirstAlgVar (A, alpha) || hasFirstAlgVar (B, alpha))
    1216   {
    1217     int degAx= degree (A, 1);
    1218     int degAy= degree (A, 2);
    1219     int degBx= degree (B, 1);
    1220     int degBy= degree (B, 2);
    1221     int d1= degAx + degBx + 1;
    1222     int d2= tmax (degAy, degBy);
    1223     zz_p::init (getCharacteristic());
    1224     zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
    1225     zz_pE::init (NTLMipo);
    1226 
    1227     int degMipo= degree (getMipo (alpha));
    1228     if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy) &&
    1229         (2*degAy > degree (M)))
    1230       return mulMod2NTLFqReci (A, B, M, alpha);
    1231 
    1232     zz_pEX NTLA= kronSub (A, d1, alpha);
    1233     zz_pEX NTLB= kronSub (B, d1, alpha);
    1234 
    1235     int k= d1*degree (M);
    1236 
    1237     MulTrunc (NTLA, NTLA, NTLB, (long) k);
    1238 
    1239     A= reverseSubst (NTLA, d1, alpha);
    1240 
    1241     return A;
    1242   }
    1243   else
    1244 #ifdef HAVE_FLINT
    1245     return mulMod2FLINTFp (A, B, M);
    1246 #else
    1247     return mulMod2NTLFp (A, B, M);
    1248 #endif
    1249 }
    1250 
    1251 #ifdef HAVE_FLINT
    1252 void
    1253 kronSubRecipro (nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm& A,
    1254                 int d)
     647kronSubReciproFp (nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm& A,
     648                  int d)
    1255649{
    1256650  int degAy= degree (A);
    … …  
    1291685
    1292686void
    1293 kronSubRecipro (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
    1294                 int d)
     687kronSubReciproQ (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
     688                 int d)
    1295689{
    1296690  int degAy= degree (A);
    … …  
    1372766
    1373767CanonicalForm
    1374 reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
     768reverseSubstReciproFp (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
    1375769{
    1376770  Variable y= Variable (2);
    … …  
    1490884
    1491885CanonicalForm
    1492 reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
     886reverseSubstReciproQ (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
    1493887{
    1494888  Variable y= Variable (2);
    … …  
    16741068  nmod_poly_init_preinv (F1, getCharacteristic(), ninv);
    16751069  nmod_poly_init_preinv (F2, getCharacteristic(), ninv);
    1676   kronSubRecipro (F1, F2, F, d1);
     1070  kronSubReciproFp (F1, F2, F, d1);
    16771071
    16781072  nmod_poly_t G1, G2;
    16791073  nmod_poly_init_preinv (G1, getCharacteristic(), ninv);
    16801074  nmod_poly_init_preinv (G2, getCharacteristic(), ninv);
    1681   kronSubRecipro (G1, G2, G, d1);
     1075  kronSubReciproFp (G1, G2, G, d1);
    16821076
    16831077  int k= d1*degree (M);
    … …  
    16961090
    16971091
    1698   CanonicalForm result= reverseSubst (F1, F2, d1, d2);
     1092  CanonicalForm result= reverseSubstReciproFp (F1, F2, d1, d2);
    16991093
    17001094  nmod_poly_clear (F1);
    … …  
    17501144  fmpz_poly_init (F1);
    17511145  fmpz_poly_init (F2);
    1752   kronSubRecipro (F1, F2, F, d1);
     1146  kronSubReciproQ (F1, F2, F, d1);
    17531147
    17541148  fmpz_poly_t G1, G2;
    17551149  fmpz_poly_init (G1);
    17561150  fmpz_poly_init (G2);
    1757   kronSubRecipro (G1, G2, G, d1);
     1151  kronSubReciproQ (G1, G2, G, d1);
    17581152
    17591153  int k= d1*degree (M);
    … …  
    17711165  int d2= tmax (fmpz_poly_degree (F2)/d1, fmpz_poly_degree (F1)/d1);
    17721166
    1773   CanonicalForm result= reverseSubst (F1, F2, d1, d2);
     1167  CanonicalForm result= reverseSubstReciproQ (F1, F2, d1, d2);
    17741168
    17751169  fmpz_poly_clear (F1);
    … …  
    18111205
    18121206#endif
     1207
     1208zz_pX kronSubFp (const CanonicalForm& A, int d)
     1209{
     1210  int degAy= degree (A);
     1211  zz_pX result;
     1212  result.rep.SetLength (d*(degAy + 1));
     1213
     1214  zz_p *resultp;
     1215  resultp= result.rep.elts();
     1216  zz_pX buf;
     1217  zz_p *bufp;
     1218  int j, k, bufRepLength;
     1219
     1220  for (CFIterator i= A; i.hasTerms(); i++)
     1221  {
     1222    if (i.coeff().inCoeffDomain())
     1223      buf= convertFacCF2NTLzzpX (i.coeff());
     1224    else
     1225      buf= convertFacCF2NTLzzpX (i.coeff());
     1226
     1227    k= i.exp()*d;
     1228    bufp= buf.rep.elts();
     1229    bufRepLength= (int) buf.rep.length();
     1230    for (j= 0; j < bufRepLength; j++)
     1231      resultp [j + k]= bufp [j];
     1232  }
     1233  result.normalize();
     1234
     1235  return result;
     1236}
     1237
     1238zz_pEX kronSubFq (const CanonicalForm& A, int d, const Variable& alpha)
     1239{
     1240  int degAy= degree (A);
     1241  zz_pEX result;
     1242  result.rep.SetLength (d*(degAy + 1));
     1243
     1244  Variable v;
     1245  zz_pE *resultp;
     1246  resultp= result.rep.elts();
     1247  zz_pEX buf1;
     1248  zz_pE *buf1p;
     1249  zz_pX buf2;
     1250  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
     1251  int j, k, buf1RepLength;
     1252
     1253  for (CFIterator i= A; i.hasTerms(); i++)
     1254  {
     1255    if (i.coeff().inCoeffDomain())
     1256    {
     1257      buf2= convertFacCF2NTLzzpX (i.coeff());
     1258      buf1= to_zz_pEX (to_zz_pE (buf2));
     1259    }
     1260    else
     1261      buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
     1262
     1263    k= i.exp()*d;
     1264    buf1p= buf1.rep.elts();
     1265    buf1RepLength= (int) buf1.rep.length();
     1266    for (j= 0; j < buf1RepLength; j++)
     1267      resultp [j + k]= buf1p [j];
     1268  }
     1269  result.normalize();
     1270
     1271  return result;
     1272}
     1273
     1274void
     1275kronSubReciproFq (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
     1276                  const Variable& alpha)
     1277{
     1278  int degAy= degree (A);
     1279  subA1.rep.SetLength ((long) d*(degAy + 2));
     1280  subA2.rep.SetLength ((long) d*(degAy + 2));
     1281
     1282  Variable v;
     1283  zz_pE *subA1p;
     1284  zz_pE *subA2p;
     1285  subA1p= subA1.rep.elts();
     1286  subA2p= subA2.rep.elts();
     1287  zz_pEX buf;
     1288  zz_pE *bufp;
     1289  zz_pX buf2;
     1290  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
     1291  int j, k, kk, bufRepLength;
     1292
     1293  for (CFIterator i= A; i.hasTerms(); i++)
     1294  {
     1295    if (i.coeff().inCoeffDomain())
     1296    {
     1297      buf2= convertFacCF2NTLzzpX (i.coeff());
     1298      buf= to_zz_pEX (to_zz_pE (buf2));
     1299    }
     1300    else
     1301      buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
     1302
     1303    k= i.exp()*d;
     1304    kk= (degAy - i.exp())*d;
     1305    bufp= buf.rep.elts();
     1306    bufRepLength= (int) buf.rep.length();
     1307    for (j= 0; j < bufRepLength; j++)
     1308    {
     1309      subA1p [j + k] += bufp [j];
     1310      subA2p [j + kk] += bufp [j];
     1311    }
     1312  }
     1313  subA1.normalize();
     1314  subA2.normalize();
     1315}
     1316
     1317void
     1318kronSubReciproFp (zz_pX& subA1, zz_pX& subA2, const CanonicalForm& A, int d)
     1319{
     1320  int degAy= degree (A);
     1321  subA1.rep.SetLength ((long) d*(degAy + 2));
     1322  subA2.rep.SetLength ((long) d*(degAy + 2));
     1323
     1324  zz_p *subA1p;
     1325  zz_p *subA2p;
     1326  subA1p= subA1.rep.elts();
     1327  subA2p= subA2.rep.elts();
     1328  zz_pX buf;
     1329  zz_p *bufp;
     1330  int j, k, kk, bufRepLength;
     1331
     1332  for (CFIterator i= A; i.hasTerms(); i++)
     1333  {
     1334    buf= convertFacCF2NTLzzpX (i.coeff());
     1335
     1336    k= i.exp()*d;
     1337    kk= (degAy - i.exp())*d;
     1338    bufp= buf.rep.elts();
     1339    bufRepLength= (int) buf.rep.length();
     1340    for (j= 0; j < bufRepLength; j++)
     1341    {
     1342      subA1p [j + k] += bufp [j];
     1343      subA2p [j + kk] += bufp [j];
     1344    }
     1345  }
     1346  subA1.normalize();
     1347  subA2.normalize();
     1348}
     1349
     1350CanonicalForm
     1351reverseSubstReciproFq (const zz_pEX& F, const zz_pEX& G, int d, int k,
     1352                       const Variable& alpha)
     1353{
     1354  Variable y= Variable (2);
     1355  Variable x= Variable (1);
     1356
     1357  zz_pEX f= F;
     1358  zz_pEX g= G;
     1359  int degf= deg(f);
     1360  int degg= deg(g);
     1361
     1362  zz_pEX buf1;
     1363  zz_pEX buf2;
     1364  zz_pEX buf3;
     1365
     1366  zz_pE *buf1p;
     1367  zz_pE *buf2p;
     1368  zz_pE *buf3p;
     1369  if (f.rep.length() < (long) d*(k+1)) //zero padding
     1370    f.rep.SetLength ((long)d*(k+1));
     1371
     1372  zz_pE *gp= g.rep.elts();
     1373  zz_pE *fp= f.rep.elts();
     1374  CanonicalForm result= 0;
     1375  int i= 0;
     1376  int lf= 0;
     1377  int lg= d*k;
     1378  int degfSubLf= degf;
     1379  int deggSubLg= degg-lg;
     1380  int repLengthBuf2, repLengthBuf1, ind, tmp;
     1381  zz_pE zzpEZero= zz_pE();
     1382
     1383  while (degf >= lf || lg >= 0)
     1384  {
     1385    if (degfSubLf >= d)
     1386      repLengthBuf1= d;
     1387    else if (degfSubLf < 0)
     1388      repLengthBuf1= 0;
     1389    else
     1390      repLengthBuf1= degfSubLf + 1;
     1391    buf1.rep.SetLength((long) repLengthBuf1);
     1392
     1393    buf1p= buf1.rep.elts();
     1394    for (ind= 0; ind < repLengthBuf1; ind++)
     1395      buf1p [ind]= fp [ind + lf];
     1396    buf1.normalize();
     1397
     1398    repLengthBuf1= buf1.rep.length();
     1399
     1400    if (deggSubLg >= d - 1)
     1401      repLengthBuf2= d - 1;
     1402    else if (deggSubLg < 0)
     1403      repLengthBuf2= 0;
     1404    else
     1405      repLengthBuf2= deggSubLg + 1;
     1406
     1407    buf2.rep.SetLength ((long) repLengthBuf2);
     1408    buf2p= buf2.rep.elts();
     1409    for (ind= 0; ind < repLengthBuf2; ind++)
     1410      buf2p [ind]= gp [ind + lg];
     1411    buf2.normalize();
     1412
     1413    repLengthBuf2= buf2.rep.length();
     1414
     1415    buf3.rep.SetLength((long) repLengthBuf2 + d);
     1416    buf3p= buf3.rep.elts();
     1417    buf2p= buf2.rep.elts();
     1418    buf1p= buf1.rep.elts();
     1419    for (ind= 0; ind < repLengthBuf1; ind++)
     1420      buf3p [ind]= buf1p [ind];
     1421    for (ind= repLengthBuf1; ind < d; ind++)
     1422      buf3p [ind]= zzpEZero;
     1423    for (ind= 0; ind < repLengthBuf2; ind++)
     1424      buf3p [ind + d]= buf2p [ind];
     1425    buf3.normalize();
     1426
     1427    result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
     1428    i++;
     1429
     1430
     1431    lf= i*d;
     1432    degfSubLf= degf - lf;
     1433
     1434    lg= d*(k-i);
     1435    deggSubLg= degg - lg;
     1436
     1437    buf1p= buf1.rep.elts();
     1438
     1439    if (lg >= 0 && deggSubLg > 0)
     1440    {
     1441      if (repLengthBuf2 > degfSubLf + 1)
     1442        degfSubLf= repLengthBuf2 - 1;
     1443      tmp= tmin (repLengthBuf1, deggSubLg + 1);
     1444      for (ind= 0; ind < tmp; ind++)
     1445        gp [ind + lg] -= buf1p [ind];
     1446    }
     1447
     1448    if (lg < 0)
     1449      break;
     1450
     1451    buf2p= buf2.rep.elts();
     1452    if (degfSubLf >= 0)
     1453    {
     1454      for (ind= 0; ind < repLengthBuf2; ind++)
     1455        fp [ind + lf] -= buf2p [ind];
     1456    }
     1457  }
     1458
     1459  return result;
     1460}
     1461
     1462CanonicalForm
     1463reverseSubstReciproFp (const zz_pX& F, const zz_pX& G, int d, int k)
     1464{
     1465  Variable y= Variable (2);
     1466  Variable x= Variable (1);
     1467
     1468  zz_pX f= F;
     1469  zz_pX g= G;
     1470  int degf= deg(f);
     1471  int degg= deg(g);
     1472
     1473  zz_pX buf1;
     1474  zz_pX buf2;
     1475  zz_pX buf3;
     1476
     1477  zz_p *buf1p;
     1478  zz_p *buf2p;
     1479  zz_p *buf3p;
     1480
     1481  if (f.rep.length() < (long) d*(k+1)) //zero padding
     1482    f.rep.SetLength ((long)d*(k+1));
     1483
     1484  zz_p *gp= g.rep.elts();
     1485  zz_p *fp= f.rep.elts();
     1486  CanonicalForm result= 0;
     1487  int i= 0;
     1488  int lf= 0;
     1489  int lg= d*k;
     1490  int degfSubLf= degf;
     1491  int deggSubLg= degg-lg;
     1492  int repLengthBuf2, repLengthBuf1, ind, tmp;
     1493  zz_p zzpZero= zz_p();
     1494  while (degf >= lf || lg >= 0)
     1495  {
     1496    if (degfSubLf >= d)
     1497      repLengthBuf1= d;
     1498    else if (degfSubLf < 0)
     1499      repLengthBuf1= 0;
     1500    else
     1501      repLengthBuf1= degfSubLf + 1;
     1502    buf1.rep.SetLength((long) repLengthBuf1);
     1503
     1504    buf1p= buf1.rep.elts();
     1505    for (ind= 0; ind < repLengthBuf1; ind++)
     1506      buf1p [ind]= fp [ind + lf];
     1507    buf1.normalize();
     1508
     1509    repLengthBuf1= buf1.rep.length();
     1510
     1511    if (deggSubLg >= d - 1)
     1512      repLengthBuf2= d - 1;
     1513    else if (deggSubLg < 0)
     1514      repLengthBuf2= 0;
     1515    else
     1516      repLengthBuf2= deggSubLg + 1;
     1517
     1518    buf2.rep.SetLength ((long) repLengthBuf2);
     1519    buf2p= buf2.rep.elts();
     1520    for (ind= 0; ind < repLengthBuf2; ind++)
     1521      buf2p [ind]= gp [ind + lg];
     1522
     1523    buf2.normalize();
     1524
     1525    repLengthBuf2= buf2.rep.length();
     1526
     1527
     1528    buf3.rep.SetLength((long) repLengthBuf2 + d);
     1529    buf3p= buf3.rep.elts();
     1530    buf2p= buf2.rep.elts();
     1531    buf1p= buf1.rep.elts();
     1532    for (ind= 0; ind < repLengthBuf1; ind++)
     1533      buf3p [ind]= buf1p [ind];
     1534    for (ind= repLengthBuf1; ind < d; ind++)
     1535      buf3p [ind]= zzpZero;
     1536    for (ind= 0; ind < repLengthBuf2; ind++)
     1537      buf3p [ind + d]= buf2p [ind];
     1538    buf3.normalize();
     1539
     1540    result += convertNTLzzpX2CF (buf3, x)*power (y, i);
     1541    i++;
     1542
     1543
     1544    lf= i*d;
     1545    degfSubLf= degf - lf;
     1546
     1547    lg= d*(k-i);
     1548    deggSubLg= degg - lg;
     1549
     1550    buf1p= buf1.rep.elts();
     1551
     1552    if (lg >= 0 && deggSubLg > 0)
     1553    {
     1554      if (repLengthBuf2 > degfSubLf + 1)
     1555        degfSubLf= repLengthBuf2 - 1;
     1556      tmp= tmin (repLengthBuf1, deggSubLg + 1);
     1557      for (ind= 0; ind < tmp; ind++)
     1558        gp [ind + lg] -= buf1p [ind];
     1559    }
     1560    if (lg < 0)
     1561      break;
     1562
     1563    buf2p= buf2.rep.elts();
     1564    if (degfSubLf >= 0)
     1565    {
     1566      for (ind= 0; ind < repLengthBuf2; ind++)
     1567        fp [ind + lf] -= buf2p [ind];
     1568    }
     1569  }
     1570
     1571  return result;
     1572}
     1573
     1574CanonicalForm reverseSubstFq (const zz_pEX& F, int d, const Variable& alpha)
     1575{
     1576  Variable y= Variable (2);
     1577  Variable x= Variable (1);
     1578
     1579  zz_pEX f= F;
     1580  zz_pE *fp= f.rep.elts();
     1581
     1582  zz_pEX buf;
     1583  zz_pE *bufp;
     1584  CanonicalForm result= 0;
     1585  int i= 0;
     1586  int degf= deg(f);
     1587  int k= 0;
     1588  int degfSubK, repLength, j;
     1589  while (degf >= k)
     1590  {
     1591    degfSubK= degf - k;
     1592    if (degfSubK >= d)
     1593      repLength= d;
     1594    else
     1595      repLength= degfSubK + 1;
     1596
     1597    buf.rep.SetLength ((long) repLength);
     1598    bufp= buf.rep.elts();
     1599    for (j= 0; j < repLength; j++)
     1600      bufp [j]= fp [j + k];
     1601    buf.normalize();
     1602
     1603    result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
     1604    i++;
     1605    k= d*i;
     1606  }
     1607
     1608  return result;
     1609}
     1610
     1611CanonicalForm reverseSubstFp (const zz_pX& F, int d)
     1612{
     1613  Variable y= Variable (2);
     1614  Variable x= Variable (1);
     1615
     1616  zz_pX f= F;
     1617  zz_p *fp= f.rep.elts();
     1618
     1619  zz_pX buf;
     1620  zz_p *bufp;
     1621  CanonicalForm result= 0;
     1622  int i= 0;
     1623  int degf= deg(f);
     1624  int k= 0;
     1625  int degfSubK, repLength, j;
     1626  while (degf >= k)
     1627  {
     1628    degfSubK= degf - k;
     1629    if (degfSubK >= d)
     1630      repLength= d;
     1631    else
     1632      repLength= degfSubK + 1;
     1633
     1634    buf.rep.SetLength ((long) repLength);
     1635    bufp= buf.rep.elts();
     1636    for (j= 0; j < repLength; j++)
     1637      bufp [j]= fp [j + k];
     1638    buf.normalize();
     1639
     1640    result += convertNTLzzpX2CF (buf, x)*power (y, i);
     1641    i++;
     1642    k= d*i;
     1643  }
     1644
     1645  return result;
     1646}
     1647
     1648// assumes input to be reduced mod M and to be an element of Fq not Fp
     1649CanonicalForm
     1650mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
     1651                  CanonicalForm& M)
     1652{
     1653  int d1= degree (F, 1) + degree (G, 1) + 1;
     1654  d1 /= 2;
     1655  d1 += 1;
     1656
     1657  zz_pX F1, F2;
     1658  kronSubReciproFp (F1, F2, F, d1);
     1659  zz_pX G1, G2;
     1660  kronSubReciproFp (G1, G2, G, d1);
     1661
     1662  int k= d1*degree (M);
     1663  MulTrunc (F1, F1, G1, (long) k);
     1664
     1665  int degtailF= degree (tailcoeff (F), 1);
     1666  int degtailG= degree (tailcoeff (G), 1);
     1667  int taildegF= taildegree (F);
     1668  int taildegG= taildegree (G);
     1669  int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
     1670
     1671  reverse (F2, F2);
     1672  reverse (G2, G2);
     1673  MulTrunc (F2, F2, G2, b + 1);
     1674  reverse (F2, F2, b);
     1675
     1676  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
     1677  return reverseSubstReciproFp (F1, F2, d1, d2);
     1678}
     1679
     1680//Kronecker substitution
     1681CanonicalForm
     1682mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
     1683              CanonicalForm& M)
     1684{
     1685  CanonicalForm A= F;
     1686  CanonicalForm B= G;
     1687
     1688  int degAx= degree (A, 1);
     1689  int degAy= degree (A, 2);
     1690  int degBx= degree (B, 1);
     1691  int degBy= degree (B, 2);
     1692  int d1= degAx + 1 + degBx;
     1693  int d2= tmax (degAy, degBy);
     1694
     1695  if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
     1696    return mulMod2NTLFpReci (A, B, M);
     1697
     1698  zz_pX NTLA= kronSubFp (A, d1);
     1699  zz_pX NTLB= kronSubFp (B, d1);
     1700
     1701  int k= d1*degree (M);
     1702  MulTrunc (NTLA, NTLA, NTLB, (long) k);
     1703
     1704  A= reverseSubstFp (NTLA, d1);
     1705
     1706  return A;
     1707}
     1708
     1709// assumes input to be reduced mod M and to be an element of Fq not Fp
     1710CanonicalForm
     1711mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
     1712                  CanonicalForm& M, const Variable& alpha)
     1713{
     1714  int d1= degree (F, 1) + degree (G, 1) + 1;
     1715  d1 /= 2;
     1716  d1 += 1;
     1717
     1718  zz_pEX F1, F2;
     1719  kronSubReciproFq (F1, F2, F, d1, alpha);
     1720  zz_pEX G1, G2;
     1721  kronSubReciproFq (G1, G2, G, d1, alpha);
     1722
     1723  int k= d1*degree (M);
     1724  MulTrunc (F1, F1, G1, (long) k);
     1725
     1726  int degtailF= degree (tailcoeff (F), 1);
     1727  int degtailG= degree (tailcoeff (G), 1);
     1728  int taildegF= taildegree (F);
     1729  int taildegG= taildegree (G);
     1730  int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
     1731
     1732  reverse (F2, F2);
     1733  reverse (G2, G2);
     1734  MulTrunc (F2, F2, G2, b + 1);
     1735  reverse (F2, F2, b);
     1736
     1737  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
     1738  return reverseSubstReciproFq (F1, F2, d1, d2, alpha);
     1739}
     1740
     1741#ifdef HAVE_FLINT
     1742CanonicalForm
     1743mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
     1744                CanonicalForm& M);
     1745#endif
     1746
     1747CanonicalForm
     1748mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
     1749              CanonicalForm& M)
     1750{
     1751  Variable alpha;
     1752  CanonicalForm A= F;
     1753  CanonicalForm B= G;
     1754
     1755  if (hasFirstAlgVar (A, alpha) || hasFirstAlgVar (B, alpha))
     1756  {
     1757    int degAx= degree (A, 1);
     1758    int degAy= degree (A, 2);
     1759    int degBx= degree (B, 1);
     1760    int degBy= degree (B, 2);
     1761    int d1= degAx + degBx + 1;
     1762    int d2= tmax (degAy, degBy);
     1763    zz_p::init (getCharacteristic());
     1764    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
     1765    zz_pE::init (NTLMipo);
     1766
     1767    int degMipo= degree (getMipo (alpha));
     1768    if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy) &&
     1769        (2*degAy > degree (M)))
     1770      return mulMod2NTLFqReci (A, B, M, alpha);
     1771
     1772    zz_pEX NTLA= kronSubFq (A, d1, alpha);
     1773    zz_pEX NTLB= kronSubFq (B, d1, alpha);
     1774
     1775    int k= d1*degree (M);
     1776
     1777    MulTrunc (NTLA, NTLA, NTLB, (long) k);
     1778
     1779    A= reverseSubstFq (NTLA, d1, alpha);
     1780
     1781    return A;
     1782  }
     1783  else
     1784#ifdef HAVE_FLINT
     1785    return mulMod2FLINTFp (A, B, M);
     1786#else
     1787    return mulMod2NTLFp (A, B, M);
     1788#endif
     1789}
    18131790
    18141791CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
    … …  
    18861863  DEBOUTLN (cerr, "fatal end in mulMod2");
    18871864}
     1865
     1866// end bivariate polys
     1867//**********************
     1868// multivariate polys
    18881869
    18891870CanonicalForm mod (const CanonicalForm& F, const CFList& M)
    … …  
    20312012}
    20322013
     2014// end multivariate polys
     2015//***************************
     2016// division
     2017
    20332018CanonicalForm reverse (const CanonicalForm& F, int d)
    20342019{
    … …  
    25102495}
    25112496
     2497// end division
     2498
    25122499#endif
Note: See TracChangeset for help on using the changeset viewer.