Changeset 0e2e23 in git


Ignore:
Timestamp:
Feb 9, 2012, 4:47:39 PM (11 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
Children:
81d96cef80a045a0f0adf2d6e851d67d9d72074e
Parents:
deb0f602cd605c4cc65110a3abe2c5eeb3b47945
git-author:
Martin Lee <martinlee84@web.de>2012-02-09 16:47:39+01:00
git-committer:
Martin Lee <martinlee84@web.de>2012-04-04 14:42:25+02:00
Message:
chg: separated multiplication and Hensel lifting functions
Location:
factory
Files:
2 added
7 edited

Legend:

Unmodified
Added
Removed

  • factory/Makefile.am

    rdeb0f60 r0e2e23  
    7676                facHensel.cc \
    7777                facIrredTest.cc \
     78                facMul.cc \
    7879                facSparseHensel.cc \
    7980                fieldGCD.cc \
    … …  
    160161                facHensel.h \
    161162                facIrredTest.h \
     163                facMul.h \
    162164                facSparseHensel.h \
    163165                fieldGCD.h \

  • factory/facBivar.cc

    rdeb0f60 r0e2e23  
    2323#include "facBivar.h"
    2424#include "facHensel.h"
     25#include "facMul.h"
    2526
    2627#ifdef HAVE_NTL

  • factory/facFqBivar.cc

    rdeb0f60 r0e2e23  
    3030#include "cf_random.h"
    3131#include "facHensel.h"
     32#include "facMul.h"
    3233#include "cf_map.h"
    3334#include "cf_gcd_smallp.h"

  • factory/facFqBivarUtil.cc

    rdeb0f60 r0e2e23  
    2727#include "cfNewtonPolygon.h"
    2828#include "facHensel.h"
     29#include "facMul.h"
    2930
    3031

  • factory/facFqFactorize.cc

    rdeb0f60 r0e2e23  
    3232#include "algext.h"
    3333#include "facSparseHensel.h"
     34#include "facMul.h"
    3435
    3536#ifdef HAVE_NTL

  • factory/facHensel.cc

    rdeb0f60 r0e2e23  
    44/** @file facHensel.cc
    55 *
    6  * This file implements functions to lift factors via Hensel lifting and
    7  * functions for modular multiplication and division with remainder.
     6 * This file implements functions to lift factors via Hensel lifting.
    87 *
    98 * ABSTRACT: Hensel lifting is described in "Efficient Multivariate
    … …  
    2726
    2827#include "facHensel.h"
     28#include "facMul.h"
    2929#include "cf_util.h"
    3030#include "fac_util.h"
    … …  
    3535#include <NTL/lzz_pEX.h>
    3636#include "NTLconvert.h"
    37 
    38 #ifdef HAVE_FLINT
    39 #include "FLINTconvert.h"
    40 #endif
    41 
    42 #ifdef HAVE_FLINT
    43 void kronSub (fmpz_poly_t result, const CanonicalForm& A, int d)
    44 {
    45   int degAy= degree (A);
    46   fmpz_poly_init2 (result, d*(degAy + 1));
    47   _fmpz_poly_set_length (result, d*(degAy + 1));
    48   CFIterator j;
    49   for (CFIterator i= A; i.hasTerms(); i++)
    50   {
    51     if (i.coeff().inBaseDomain())
    52       convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d), i.coeff());
    53     else
    54       for (j= i.coeff(); j.hasTerms(); j++)
    55         convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d+j.exp()),
    56                         j.coeff());
    57   }
    58   _fmpz_poly_normalise(result);
    59 }
    60 
    61 
    62 CanonicalForm
    63 reverseSubst (const fmpz_poly_t F, int d, const Variable& alpha,
    64               const CanonicalForm& den)
    65 {
    66   Variable x= Variable (1);
    67 
    68   CanonicalForm result= 0;
    69   int i= 0;
    70   int degf= fmpz_poly_degree (F);
    71   int k= 0;
    72   int degfSubK;
    73   int repLength, j;
    74   CanonicalForm coeff;
    75   fmpz* tmp;
    76   while (degf >= k)
    77   {
    78     coeff= 0;
    79     degfSubK= degf - k;
    80     if (degfSubK >= d)
    81       repLength= d;
    82     else
    83       repLength= degfSubK + 1;
    84 
    85     for (j= 0; j < repLength; j++)
    86     {
    87       tmp= fmpz_poly_get_coeff_ptr (F, j+k);
    88       if (!fmpz_is_zero (tmp))
    89       {
    90         CanonicalForm ff= convertFmpz2CF (tmp)/den;
    91         coeff += ff*power (alpha, j);
    92       }
    93     }
    94     result += coeff*power (x, i);
    95     i++;
    96     k= d*i;
    97   }
    98   return result;
    99 }
    100 
    101 CanonicalForm
    102 mulFLINTQa (const CanonicalForm& F, const CanonicalForm& G,
    103             const Variable& alpha)
    104 {
    105   CanonicalForm A= F;
    106   CanonicalForm B= G;
    107 
    108   CanonicalForm denA= bCommonDen (A);
    109   CanonicalForm denB= bCommonDen (B);
    110 
    111   A *= denA;
    112   B *= denB;
    113   int degAa= degree (A, alpha);
    114   int degBa= degree (B, alpha);
    115   int d= degAa + 1 + degBa;
    116 
    117   fmpz_poly_t FLINTA,FLINTB;
    118   fmpz_poly_init (FLINTA);
    119   fmpz_poly_init (FLINTB);
    120   kronSub (FLINTA, A, d);
    121   kronSub (FLINTB, B, d);
    122 
    123   fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
    124 
    125   denA *= denB;
    126   A= reverseSubst (FLINTA, d, alpha, denA);
    127 
    128   fmpz_poly_clear (FLINTA);
    129   fmpz_poly_clear (FLINTB);
    130   return A;
    131 }
    132 
    133 CanonicalForm
    134 mulFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
    135 {
    136   CanonicalForm A= F;
    137   CanonicalForm B= G;
    138 
    139   CanonicalForm denA= bCommonDen (A);
    140   CanonicalForm denB= bCommonDen (B);
    141 
    142   A *= denA;
    143   B *= denB;
    144   fmpz_poly_t FLINTA,FLINTB;
    145   convertFacCF2Fmpz_poly_t (FLINTA, A);
    146   convertFacCF2Fmpz_poly_t (FLINTB, B);
    147   fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
    148   denA *= denB;
    149   A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
    150   A /= denA;
    151   fmpz_poly_clear (FLINTA);
    152   fmpz_poly_clear (FLINTB);
    153 
    154   return A;
    155 }
    156 
    157 /*CanonicalForm
    158 mulFLINTQ2 (const CanonicalForm& F, const CanonicalForm& G)
    159 {
    160   CanonicalForm A= F;
    161   CanonicalForm B= G;
    162 
    163   fmpq_poly_t FLINTA,FLINTB;
    164   convertFacCF2Fmpq_poly_t (FLINTA, A);
    165   convertFacCF2Fmpq_poly_t (FLINTB, B);
    166 
    167   fmpq_poly_mul (FLINTA, FLINTA, FLINTB);
    168   A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
    169   fmpq_poly_clear (FLINTA);
    170   fmpq_poly_clear (FLINTB);
    171   return A;
    172 }*/
    173 
    174 CanonicalForm
    175 divFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
    176 {
    177   CanonicalForm A= F;
    178   CanonicalForm B= G;
    179 
    180   fmpq_poly_t FLINTA,FLINTB;
    181   convertFacCF2Fmpq_poly_t (FLINTA, A);
    182   convertFacCF2Fmpq_poly_t (FLINTB, B);
    183 
    184   fmpq_poly_div (FLINTA, FLINTA, FLINTB);
    185   A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
    186 
    187   fmpq_poly_clear (FLINTA);
    188   fmpq_poly_clear (FLINTB);
    189   return A;
    190 }
    191 
    192 CanonicalForm
    193 modFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
    194 {
    195   CanonicalForm A= F;
    196   CanonicalForm B= G;
    197 
    198   fmpq_poly_t FLINTA,FLINTB;
    199   convertFacCF2Fmpq_poly_t (FLINTA, A);
    200   convertFacCF2Fmpq_poly_t (FLINTB, B);
    201 
    202   fmpq_poly_rem (FLINTA, FLINTA, FLINTB);
    203   A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
    204 
    205   fmpq_poly_clear (FLINTA);
    206   fmpq_poly_clear (FLINTB);
    207   return A;
    208 }
    209 
    210 CanonicalForm
    211 mulFLINTQaTrunc (const CanonicalForm& F, const CanonicalForm& G,
    212                  const Variable& alpha, int m)
    213 {
    214   CanonicalForm A= F;
    215   CanonicalForm B= G;
    216 
    217   CanonicalForm denA= bCommonDen (A);
    218   CanonicalForm denB= bCommonDen (B);
    219 
    220   A *= denA;
    221   B *= denB;
    222 
    223   int degAa= degree (A, alpha);
    224   int degBa= degree (B, alpha);
    225   int d= degAa + 1 + degBa;
    226 
    227   fmpz_poly_t FLINTA,FLINTB;
    228   fmpz_poly_init (FLINTA);
    229   fmpz_poly_init (FLINTB);
    230   kronSub (FLINTA, A, d);
    231   kronSub (FLINTB, B, d);
    232 
    233   int k= d*m;
    234   fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, k);
    235 
    236   denA *= denB;
    237   A= reverseSubst (FLINTA, d, alpha, denA);
    238   fmpz_poly_clear (FLINTA);
    239   fmpz_poly_clear (FLINTB);
    240   return A;
    241 }
    242 
    243 CanonicalForm
    244 mulFLINTQTrunc (const CanonicalForm& F, const CanonicalForm& G, int m)
    245 {
    246   if (F.inCoeffDomain() || G.inCoeffDomain())
    247     return mod (F*G, power (Variable (1), m));
    248   Variable alpha;
    249   if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
    250     return mulFLINTQaTrunc (F, G, alpha, m);
    251 
    252   CanonicalForm A= F;
    253   CanonicalForm B= G;
    254 
    255   CanonicalForm denA= bCommonDen (A);
    256   CanonicalForm denB= bCommonDen (B);
    257 
    258   A *= denA;
    259   B *= denB;
    260   fmpz_poly_t FLINTA,FLINTB;
    261   convertFacCF2Fmpz_poly_t (FLINTA, A);
    262   convertFacCF2Fmpz_poly_t (FLINTB, B);
    263   fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, m);
    264   denA *= denB;
    265   A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
    266   A /= denA;
    267   fmpz_poly_clear (FLINTA);
    268   fmpz_poly_clear (FLINTB);
    269 
    270   return A;
    271 }
    272 
    273 CanonicalForm uniReverse (const CanonicalForm& F, int d)
    274 {
    275   if (d == 0)
    276     return F;
    277   if (F.inCoeffDomain())
    278     return F*power (Variable (1),d);
    279   Variable x= Variable (1);
    280   CanonicalForm result= 0;
    281   CFIterator i= F;
    282   while (d - i.exp() < 0)
    283     i++;
    284 
    285   for (; i.hasTerms() && (d - i.exp() >= 0); i++)
    286     result += i.coeff()*power (x, d - i.exp());
    287   return result;
    288 }
    289 
    290 CanonicalForm
    291 newtonInverse (const CanonicalForm& F, const int n)
    292 {
    293   int l= ilog2(n);
    294 
    295   CanonicalForm g= F [0];
    296 
    297   ASSERT (!g.isZero(), "expected a unit");
    298 
    299   if (!g.isOne())
    300     g = 1/g;
    301   Variable x= Variable (1);
    302   CanonicalForm result;
    303   int exp= 0;
    304   if (n & 1)
    305   {
    306     result= g;
    307     exp= 1;
    308   }
    309   CanonicalForm h;
    310 
    311   for (int i= 1; i <= l; i++)
    312   {
    313     h= mulNTL (g, mod (F, power (x, (1 << i))));
    314     h= mod (h, power (x, (1 << i)) - 1);
    315     h= div (h, power (x, (1 << (i - 1))));
    316     g -= power (x, (1 << (i - 1)))*
    317          mulFLINTQTrunc (g, h, 1 << (i-1));
    318 
    319     if (n & (1 << i))
    320     {
    321       if (exp)
    322       {
    323         h= mulNTL (result, mod (F, power (x, exp + (1 << i))));
    324         h= mod (h, power (x, exp + (1 << i)) - 1);
    325         h= div (h, power (x, exp));
    326         result -= power(x, exp)*mulFLINTQTrunc (g, h, 1 << i);
    327         exp += (1 << i);
    328       }
    329       else
    330       {
    331         exp= (1 << i);
    332         result= g;
    333       }
    334     }
    335   }
    336 
    337   return result;
    338 }
    339 
    340 void
    341 newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
    342               CanonicalForm& R)
    343 {
    344   CanonicalForm A= F;
    345   CanonicalForm B= G;
    346   Variable x= Variable (1);
    347   int degA= degree (A, x);
    348   int degB= degree (B, x);
    349   int m= degA - degB;
    350 
    351   if (m < 0)
    352   {
    353     R= A;
    354     Q= 0;
    355     return;
    356   }
    357 
    358   if (degB <= 1)
    359     divrem (A, B, Q, R);
    360   else
    361   {
    362     R= uniReverse (A, degA);
    363 
    364     CanonicalForm revB= uniReverse (B, degB);
    365     CanonicalForm buf= revB;
    366     revB= newtonInverse (revB, m + 1);
    367     Q= mulFLINTQTrunc (R, revB, m + 1);
    368     Q= uniReverse (Q, m);
    369 
    370     R= A - mulNTL (Q, B);
    371   }
    372 }
    373 
    374 void
    375 newtonDiv (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q)
    376 {
    377   CanonicalForm A= F;
    378   CanonicalForm B= G;
    379   Variable x= Variable (1);
    380   int degA= degree (A, x);
    381   int degB= degree (B, x);
    382   int m= degA - degB;
    383 
    384   if (m < 0)
    385   {
    386     Q= 0;
    387     return;
    388   }
    389 
    390   if (degB <= 1)
    391     Q= div (A, B);
    392   else
    393   {
    394     CanonicalForm R= uniReverse (A, degA);
    395 
    396     CanonicalForm revB= uniReverse (B, degB);
    397     revB= newtonInverse (revB, m + 1);
    398     Q= mulFLINTQTrunc (R, revB, m + 1);
    399     Q= uniReverse (Q, m);
    400   }
    401 }
    402 
    403 #endif
    40437
    40538static
    … …  
    697330  if (save_rat) Off(SW_RATIONAL);
    698331  return result;
    699 }
    700 
    701 CanonicalForm
    702 mulNTL (const CanonicalForm& F, const CanonicalForm& G)
    703 {
    704   if (F.inCoeffDomain() || G.inCoeffDomain() || getCharacteristic() == 0)
    705   {
    706     Variable alpha;
    707 #ifdef HAVE_FLINT
    708     if ((!F.inCoeffDomain() && !G.inCoeffDomain()) &&
    709         (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha)))
    710     {
    711       CanonicalForm result= mulFLINTQa (F, G, alpha);
    712       return result;
    713     }
    714     else if (!F.inCoeffDomain() && !G.inCoeffDomain())
    715       return mulFLINTQ (F, G);
    716 #endif
    717     return F*G;
    718   }
    719   ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
    720   ASSERT (F.level() == G.level(), "expected polys of same level");
    721   if (CFFactory::gettype() == GaloisFieldDomain)
    722     return F*G;
    723   zz_p::init (getCharacteristic());
    724   Variable alpha;
    725   CanonicalForm result;
    726   if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
    727   {
    728     zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
    729     zz_pE::init (NTLMipo);
    730     zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
    731     zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
    732     mul (NTLF, NTLF, NTLG);
    733     result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
    734   }
    735   else
    736   {
    737 #ifdef HAVE_FLINT
    738     nmod_poly_t FLINTF, FLINTG;
    739     convertFacCF2nmod_poly_t (FLINTF, F);
    740     convertFacCF2nmod_poly_t (FLINTG, G);
    741     nmod_poly_mul (FLINTF, FLINTF, FLINTG);
    742     result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
    743     nmod_poly_clear (FLINTF);
    744     nmod_poly_clear (FLINTG);
    745 #else
    746     zz_pX NTLF= convertFacCF2NTLzzpX (F);
    747     zz_pX NTLG= convertFacCF2NTLzzpX (G);
    748     mul (NTLF, NTLF, NTLG);
    749     result= convertNTLzzpX2CF(NTLF, F.mvar());
    750 #endif
    751   }
    752   return result;
    753 }
    754 
    755 CanonicalForm
    756 modNTL (const CanonicalForm& F, const CanonicalForm& G)
    757 {
    758   if (F.inCoeffDomain() && G.isUnivariate())
    759     return F;
    760   else if (F.inCoeffDomain() && G.inCoeffDomain())
    761     return mod (F, G);
    762   else if (F.isUnivariate() && G.inCoeffDomain())
    763     return mod (F,G);
    764 
    765   if (getCharacteristic() == 0)
    766   {
    767 #ifdef HAVE_FLINT
    768     Variable alpha;
    769     if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
    770       return modFLINTQ (F, G);
    771     else
    772     {
    773       CanonicalForm Q, R;
    774       newtonDivrem (F, G, Q, R);
    775       return R;
    776     }
    777 #else
    778     return mod (F, G);
    779 #endif
    780   }
    781 
    782   ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
    783   ASSERT (F.level() == G.level(), "expected polys of same level");
    784   if (CFFactory::gettype() == GaloisFieldDomain)
    785     return mod (F, G);
    786   zz_p::init (getCharacteristic());
    787   Variable alpha;
    788   CanonicalForm result;
    789   if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
    790   {
    791     zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
    792     zz_pE::init (NTLMipo);
    793     zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
    794     zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
    795     rem (NTLF, NTLF, NTLG);
    796     result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
    797   }
    798   else
    799   {
    800 #ifdef HAVE_FLINT
    801     nmod_poly_t FLINTF, FLINTG;
    802     convertFacCF2nmod_poly_t (FLINTF, F);
    803     convertFacCF2nmod_poly_t (FLINTG, G);
    804     nmod_poly_divrem (FLINTG, FLINTF, FLINTF, FLINTG);
    805     result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
    806     nmod_poly_clear (FLINTF);
    807     nmod_poly_clear (FLINTG);
    808 #else
    809     zz_pX NTLF= convertFacCF2NTLzzpX (F);
    810     zz_pX NTLG= convertFacCF2NTLzzpX (G);
    811     rem (NTLF, NTLF, NTLG);
    812     result= convertNTLzzpX2CF(NTLF, F.mvar());
    813 #endif
    814   }
    815   return result;
    816 }
    817 
    818 CanonicalForm
    819 divNTL (const CanonicalForm& F, const CanonicalForm& G)
    820 {
    821   if (F.inCoeffDomain() && G.isUnivariate())
    822     return F;
    823   else if (F.inCoeffDomain() && G.inCoeffDomain())
    824     return div (F, G);
    825   else if (F.isUnivariate() && G.inCoeffDomain())
    826     return div (F,G);
    827 
    828   if (getCharacteristic() == 0)
    829   {
    830 #ifdef HAVE_FLINT
    831     Variable alpha;
    832     if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
    833       return divFLINTQ (F,G);
    834     else
    835     {
    836       CanonicalForm Q;
    837       newtonDiv (F, G, Q);
    838       return Q;
    839     }
    840 #else
    841     return div (F, G);
    842 #endif
    843   }
    844 
    845   ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
    846   ASSERT (F.level() == G.level(), "expected polys of same level");
    847   if (CFFactory::gettype() == GaloisFieldDomain)
    848     return div (F, G);
    849   zz_p::init (getCharacteristic());
    850   Variable alpha;
    851   CanonicalForm result;
    852   if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
    853   {
    854     zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
    855     zz_pE::init (NTLMipo);
    856     zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
    857     zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
    858     div (NTLF, NTLF, NTLG);
    859     result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
    860   }
    861   else
    862   {
    863 #ifdef HAVE_FLINT
    864     nmod_poly_t FLINTF, FLINTG;
    865     convertFacCF2nmod_poly_t (FLINTF, F);
    866     convertFacCF2nmod_poly_t (FLINTG, G);
    867     nmod_poly_div (FLINTF, FLINTF, FLINTG);
    868     result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
    869     nmod_poly_clear (FLINTF);
    870     nmod_poly_clear (FLINTG);
    871 #else
    872     zz_pX NTLF= convertFacCF2NTLzzpX (F);
    873     zz_pX NTLG= convertFacCF2NTLzzpX (G);
    874     div (NTLF, NTLF, NTLG);
    875     result= convertNTLzzpX2CF(NTLF, F.mvar());
    876 #endif
    877   }
    878   return result;
    879 }
    880 
    881 /*
    882 void
    883 divremNTL (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
    884            CanonicalForm& R)
    885 {
    886   if (F.inCoeffDomain() && G.isUnivariate())
    887   {
    888     R= F;
    889     Q= 0;
    890   }
    891   else if (F.inCoeffDomain() && G.inCoeffDomain())
    892   {
    893     divrem (F, G, Q, R);
    894     return;
    895   }
    896   else if (F.isUnivariate() && G.inCoeffDomain())
    897   {
    898     divrem (F, G, Q, R);
    899     return;
    900   }
    901 
    902   ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
    903   ASSERT (F.level() == G.level(), "expected polys of same level");
    904   if (CFFactory::gettype() == GaloisFieldDomain)
    905   {
    906     divrem (F, G, Q, R);
    907     return;
    908   }
    909   zz_p::init (getCharacteristic());
    910   Variable alpha;
    911   CanonicalForm result;
    912   if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
    913   {
    914     zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
    915     zz_pE::init (NTLMipo);
    916     zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
    917     zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
    918     zz_pEX NTLQ;
    919     zz_pEX NTLR;
    920     DivRem (NTLQ, NTLR, NTLF, NTLG);
    921     Q= convertNTLzz_pEX2CF(NTLQ, F.mvar(), alpha);
    922     R= convertNTLzz_pEX2CF(NTLR, F.mvar(), alpha);
    923     return;
    924   }
    925   else
    926   {
    927     zz_pX NTLF= convertFacCF2NTLzzpX (F);
    928     zz_pX NTLG= convertFacCF2NTLzzpX (G);
    929     zz_pX NTLQ;
    930     zz_pX NTLR;
    931     DivRem (NTLQ, NTLR, NTLF, NTLG);
    932     Q= convertNTLzzpX2CF(NTLQ, F.mvar());
    933     R= convertNTLzzpX2CF(NTLR, F.mvar());
    934     return;
    935   }
    936 }*/
    937 
    938 CanonicalForm mod (const CanonicalForm& F, const CFList& M)
    939 {
    940   CanonicalForm A= F;
    941   for (CFListIterator i= M; i.hasItem(); i++)
    942     A= mod (A, i.getItem());
    943   return A;
    944 }
    945 
    946 #ifdef HAVE_FLINT
    947 void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
    948 {
    949   int degAy= degree (A);
    950   nmod_poly_init2 (result, getCharacteristic(), d*(degAy + 1));
    951 
    952   nmod_poly_t buf;
    953 
    954   int j, k, bufRepLength;
    955   for (CFIterator i= A; i.hasTerms(); i++)
    956   {
    957     convertFacCF2nmod_poly_t (buf, i.coeff());
    958 
    959     k= i.exp()*d;
    960     bufRepLength= (int) nmod_poly_length (buf);
    961     for (j= 0; j < bufRepLength; j++)
    962       nmod_poly_set_coeff_ui (result, j + k, nmod_poly_get_coeff_ui (buf, j));
    963     nmod_poly_clear (buf);
    964   }
    965   _nmod_poly_normalise (result);
    966 }
    967 
    968 /*void kronSubQ (fmpz_poly_t result, const CanonicalForm& A, int d)
    969 {
    970   int degAy= degree (A);
    971   fmpz_poly_init2 (result, d*(degAy + 1));
    972   _fmpz_poly_set_length (result, d*(degAy+1));
    973 
    974   CFIterator j;
    975   for (CFIterator i= A; i.hasTerms(); i++)
    976   {
    977     if (i.coeff().inBas
    978     convertFacCF2Fmpz_poly_t (buf, i.coeff());
    979 
    980     int k= i.exp()*d;
    981     int bufRepLength= (int) fmpz_poly_length (buf);
    982     for (int j= 0; j < bufRepLength; j++)
    983     {
    984       fmpz_poly_get_coeff_fmpz (coeff, buf, j);
    985       fmpz_poly_set_coeff_fmpz (result, j + k, coeff);
    986     }
    987     fmpz_poly_clear (buf);
    988   }
    989   fmpz_clear (coeff);
    990   _fmpz_poly_normalise (result);
    991 }*/
    992 
    993 // A is a bivariate poly over Qa!!!!
    994 // d2= 2*deg_alpha + 1
    995 // d1= 2*deg_x*d2+1
    996 void kronSubQa (fmpq_poly_t result, const CanonicalForm& A, int d1, int d2)
    997 {
    998   int degAy= degree (A);
    999   fmpq_poly_init2 (result, d1*(degAy + 1));
    1000 
    1001   fmpq_poly_t buf;
    1002   fmpq_t coeff;
    1003 
    1004   int k, l, bufRepLength;
    1005   CFIterator j;
    1006   for (CFIterator i= A; i.hasTerms(); i++)
    1007   {
    1008     if (i.coeff().inCoeffDomain())
    1009     {
    1010       k= d1*i.exp();
    1011       convertFacCF2Fmpq_poly_t (buf, i.coeff());
    1012       bufRepLength= (int) fmpq_poly_length(buf);
    1013       for (l= 0; l < bufRepLength; l++)
    1014       {
    1015         fmpq_poly_get_coeff_fmpq (coeff, buf, l);
    1016         fmpq_poly_set_coeff_fmpq (result, l + k, coeff);
    1017       }
    1018       fmpq_poly_clear (buf);
    1019     }
    1020     else
    1021     {
    1022       for (j= i.coeff(); j.hasTerms(); j++)
    1023       {
    1024         k= d1*i.exp();
    1025         k += d2*j.exp();
    1026         convertFacCF2Fmpq_poly_t (buf, j.coeff());
    1027         bufRepLength= (int) fmpq_poly_length(buf);
    1028         for (l= 0; l < bufRepLength; l++)
    1029         {
    1030           fmpq_poly_get_coeff_fmpq (coeff, buf, l);
    1031           fmpq_poly_set_coeff_fmpq (result, k + l, coeff);
    1032         }
    1033         fmpq_poly_clear (buf);
    1034       }
    1035     }
    1036   }
    1037   fmpq_clear (coeff);
    1038   _fmpq_poly_normalise (result);
    1039 }
    1040 #endif
    1041 
    1042 zz_pX kronSubFp (const CanonicalForm& A, int d)
    1043 {
    1044   int degAy= degree (A);
    1045   zz_pX result;
    1046   result.rep.SetLength (d*(degAy + 1));
    1047 
    1048   zz_p *resultp;
    1049   resultp= result.rep.elts();
    1050   zz_pX buf;
    1051   zz_p *bufp;
    1052   int j, k, bufRepLength;
    1053 
    1054   for (CFIterator i= A; i.hasTerms(); i++)
    1055   {
    1056     if (i.coeff().inCoeffDomain())
    1057       buf= convertFacCF2NTLzzpX (i.coeff());
    1058     else
    1059       buf= convertFacCF2NTLzzpX (i.coeff());
    1060 
    1061     k= i.exp()*d;
    1062     bufp= buf.rep.elts();
    1063     bufRepLength= (int) buf.rep.length();
    1064     for (j= 0; j < bufRepLength; j++)
    1065       resultp [j + k]= bufp [j];
    1066   }
    1067   result.normalize();
    1068 
    1069   return result;
    1070 }
    1071 
    1072 zz_pEX kronSub (const CanonicalForm& A, int d, const Variable& alpha)
    1073 {
    1074   int degAy= degree (A);
    1075   zz_pEX result;
    1076   result.rep.SetLength (d*(degAy + 1));
    1077 
    1078   Variable v;
    1079   zz_pE *resultp;
    1080   resultp= result.rep.elts();
    1081   zz_pEX buf1;
    1082   zz_pE *buf1p;
    1083   zz_pX buf2;
    1084   zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
    1085   int j, k, buf1RepLength;
    1086 
    1087   for (CFIterator i= A; i.hasTerms(); i++)
    1088   {
    1089     if (i.coeff().inCoeffDomain())
    1090     {
    1091       buf2= convertFacCF2NTLzzpX (i.coeff());
    1092       buf1= to_zz_pEX (to_zz_pE (buf2));
    1093     }
    1094     else
    1095       buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
    1096 
    1097     k= i.exp()*d;
    1098     buf1p= buf1.rep.elts();
    1099     buf1RepLength= (int) buf1.rep.length();
    1100     for (j= 0; j < buf1RepLength; j++)
    1101       resultp [j + k]= buf1p [j];
    1102   }
    1103   result.normalize();
    1104 
    1105   return result;
    1106 }
    1107 
    1108 void
    1109 kronSubRecipro (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
    1110                 const Variable& alpha)
    1111 {
    1112   int degAy= degree (A);
    1113   subA1.rep.SetLength ((long) d*(degAy + 2));
    1114   subA2.rep.SetLength ((long) d*(degAy + 2));
    1115 
    1116   Variable v;
    1117   zz_pE *subA1p;
    1118   zz_pE *subA2p;
    1119   subA1p= subA1.rep.elts();
    1120   subA2p= subA2.rep.elts();
    1121   zz_pEX buf;
    1122   zz_pE *bufp;
    1123   zz_pX buf2;
    1124   zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
    1125   int j, k, kk, bufRepLength;
    1126 
    1127   for (CFIterator i= A; i.hasTerms(); i++)
    1128   {
    1129     if (i.coeff().inCoeffDomain())
    1130     {
    1131       buf2= convertFacCF2NTLzzpX (i.coeff());
    1132       buf= to_zz_pEX (to_zz_pE (buf2));
    1133     }
    1134     else
    1135       buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
    1136 
    1137     k= i.exp()*d;
    1138     kk= (degAy - i.exp())*d;
    1139     bufp= buf.rep.elts();
    1140     bufRepLength= (int) buf.rep.length();
    1141     for (j= 0; j < bufRepLength; j++)
    1142     {
    1143       subA1p [j + k] += bufp [j];
    1144       subA2p [j + kk] += bufp [j];
    1145     }
    1146   }
    1147   subA1.normalize();
    1148   subA2.normalize();
    1149 }
    1150 
    1151 void
    1152 kronSubRecipro (zz_pX& subA1, zz_pX& subA2,const CanonicalForm& A, int d)
    1153 {
    1154   int degAy= degree (A);
    1155   subA1.rep.SetLength ((long) d*(degAy + 2));
    1156   subA2.rep.SetLength ((long) d*(degAy + 2));
    1157 
    1158   zz_p *subA1p;
    1159   zz_p *subA2p;
    1160   subA1p= subA1.rep.elts();
    1161   subA2p= subA2.rep.elts();
    1162   zz_pX buf;
    1163   zz_p *bufp;
    1164   int j, k, kk, bufRepLength;
    1165 
    1166   for (CFIterator i= A; i.hasTerms(); i++)
    1167   {
    1168     buf= convertFacCF2NTLzzpX (i.coeff());
    1169 
    1170     k= i.exp()*d;
    1171     kk= (degAy - i.exp())*d;
    1172     bufp= buf.rep.elts();
    1173     bufRepLength= (int) buf.rep.length();
    1174     for (j= 0; j < bufRepLength; j++)
    1175     {
    1176       subA1p [j + k] += bufp [j];
    1177       subA2p [j + kk] += bufp [j];
    1178     }
    1179   }
    1180   subA1.normalize();
    1181   subA2.normalize();
    1182 }
    1183 
    1184 CanonicalForm
    1185 reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
    1186               const Variable& alpha)
    1187 {
    1188   Variable y= Variable (2);
    1189   Variable x= Variable (1);
    1190 
    1191   zz_pEX f= F;
    1192   zz_pEX g= G;
    1193   int degf= deg(f);
    1194   int degg= deg(g);
    1195 
    1196   zz_pEX buf1;
    1197   zz_pEX buf2;
    1198   zz_pEX buf3;
    1199 
    1200   zz_pE *buf1p;
    1201   zz_pE *buf2p;
    1202   zz_pE *buf3p;
    1203   if (f.rep.length() < (long) d*(k+1)) //zero padding
    1204     f.rep.SetLength ((long)d*(k+1));
    1205 
    1206   zz_pE *gp= g.rep.elts();
    1207   zz_pE *fp= f.rep.elts();
    1208   CanonicalForm result= 0;
    1209   int i= 0;
    1210   int lf= 0;
    1211   int lg= d*k;
    1212   int degfSubLf= degf;
    1213   int deggSubLg= degg-lg;
    1214   int repLengthBuf2, repLengthBuf1, ind, tmp;
    1215   zz_pE zzpEZero= zz_pE();
    1216 
    1217   while (degf >= lf || lg >= 0)
    1218   {
    1219     if (degfSubLf >= d)
    1220       repLengthBuf1= d;
    1221     else if (degfSubLf < 0)
    1222       repLengthBuf1= 0;
    1223     else
    1224       repLengthBuf1= degfSubLf + 1;
    1225     buf1.rep.SetLength((long) repLengthBuf1);
    1226 
    1227     buf1p= buf1.rep.elts();
    1228     for (ind= 0; ind < repLengthBuf1; ind++)
    1229       buf1p [ind]= fp [ind + lf];
    1230     buf1.normalize();
    1231 
    1232     repLengthBuf1= buf1.rep.length();
    1233 
    1234     if (deggSubLg >= d - 1)
    1235       repLengthBuf2= d - 1;
    1236     else if (deggSubLg < 0)
    1237       repLengthBuf2= 0;
    1238     else
    1239       repLengthBuf2= deggSubLg + 1;
    1240 
    1241     buf2.rep.SetLength ((long) repLengthBuf2);
    1242     buf2p= buf2.rep.elts();
    1243     for (ind= 0; ind < repLengthBuf2; ind++)
    1244       buf2p [ind]= gp [ind + lg];
    1245     buf2.normalize();
    1246 
    1247     repLengthBuf2= buf2.rep.length();
    1248 
    1249     buf3.rep.SetLength((long) repLengthBuf2 + d);
    1250     buf3p= buf3.rep.elts();
    1251     buf2p= buf2.rep.elts();
    1252     buf1p= buf1.rep.elts();
    1253     for (ind= 0; ind < repLengthBuf1; ind++)
    1254       buf3p [ind]= buf1p [ind];
    1255     for (ind= repLengthBuf1; ind < d; ind++)
    1256       buf3p [ind]= zzpEZero;
    1257     for (ind= 0; ind < repLengthBuf2; ind++)
    1258       buf3p [ind + d]= buf2p [ind];
    1259     buf3.normalize();
    1260 
    1261     result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
    1262     i++;
    1263 
    1264 
    1265     lf= i*d;
    1266     degfSubLf= degf - lf;
    1267 
    1268     lg= d*(k-i);
    1269     deggSubLg= degg - lg;
    1270 
    1271     buf1p= buf1.rep.elts();
    1272 
    1273     if (lg >= 0 && deggSubLg > 0)
    1274     {
    1275       if (repLengthBuf2 > degfSubLf + 1)
    1276         degfSubLf= repLengthBuf2 - 1;
    1277       tmp= tmin (repLengthBuf1, deggSubLg + 1);
    1278       for (ind= 0; ind < tmp; ind++)
    1279         gp [ind + lg] -= buf1p [ind];
    1280     }
    1281 
    1282     if (lg < 0)
    1283       break;
    1284 
    1285     buf2p= buf2.rep.elts();
    1286     if (degfSubLf >= 0)
    1287     {
    1288       for (ind= 0; ind < repLengthBuf2; ind++)
    1289         fp [ind + lf] -= buf2p [ind];
    1290     }
    1291   }
    1292 
    1293   return result;
    1294 }
    1295 
    1296 CanonicalForm
    1297 reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
    1298 {
    1299   Variable y= Variable (2);
    1300   Variable x= Variable (1);
    1301 
    1302   zz_pX f= F;
    1303   zz_pX g= G;
    1304   int degf= deg(f);
    1305   int degg= deg(g);
    1306 
    1307   zz_pX buf1;
    1308   zz_pX buf2;
    1309   zz_pX buf3;
    1310 
    1311   zz_p *buf1p;
    1312   zz_p *buf2p;
    1313   zz_p *buf3p;
    1314 
    1315   if (f.rep.length() < (long) d*(k+1)) //zero padding
    1316     f.rep.SetLength ((long)d*(k+1));
    1317 
    1318   zz_p *gp= g.rep.elts();
    1319   zz_p *fp= f.rep.elts();
    1320   CanonicalForm result= 0;
    1321   int i= 0;
    1322   int lf= 0;
    1323   int lg= d*k;
    1324   int degfSubLf= degf;
    1325   int deggSubLg= degg-lg;
    1326   int repLengthBuf2, repLengthBuf1, ind, tmp;
    1327   zz_p zzpZero= zz_p();
    1328   while (degf >= lf || lg >= 0)
    1329   {
    1330     if (degfSubLf >= d)
    1331       repLengthBuf1= d;
    1332     else if (degfSubLf < 0)
    1333       repLengthBuf1= 0;
    1334     else
    1335       repLengthBuf1= degfSubLf + 1;
    1336     buf1.rep.SetLength((long) repLengthBuf1);
    1337 
    1338     buf1p= buf1.rep.elts();
    1339     for (ind= 0; ind < repLengthBuf1; ind++)
    1340       buf1p [ind]= fp [ind + lf];
    1341     buf1.normalize();
    1342 
    1343     repLengthBuf1= buf1.rep.length();
    1344 
    1345     if (deggSubLg >= d - 1)
    1346       repLengthBuf2= d - 1;
    1347     else if (deggSubLg < 0)
    1348       repLengthBuf2= 0;
    1349     else
    1350       repLengthBuf2= deggSubLg + 1;
    1351 
    1352     buf2.rep.SetLength ((long) repLengthBuf2);
    1353     buf2p= buf2.rep.elts();
    1354     for (ind= 0; ind < repLengthBuf2; ind++)
    1355       buf2p [ind]= gp [ind + lg];
    1356 
    1357     buf2.normalize();
    1358 
    1359     repLengthBuf2= buf2.rep.length();
    1360 
    1361 
    1362     buf3.rep.SetLength((long) repLengthBuf2 + d);
    1363     buf3p= buf3.rep.elts();
    1364     buf2p= buf2.rep.elts();
    1365     buf1p= buf1.rep.elts();
    1366     for (ind= 0; ind < repLengthBuf1; ind++)
    1367       buf3p [ind]= buf1p [ind];
    1368     for (ind= repLengthBuf1; ind < d; ind++)
    1369       buf3p [ind]= zzpZero;
    1370     for (ind= 0; ind < repLengthBuf2; ind++)
    1371       buf3p [ind + d]= buf2p [ind];
    1372     buf3.normalize();
    1373 
    1374     result += convertNTLzzpX2CF (buf3, x)*power (y, i);
    1375     i++;
    1376 
    1377 
    1378     lf= i*d;
    1379     degfSubLf= degf - lf;
    1380 
    1381     lg= d*(k-i);
    1382     deggSubLg= degg - lg;
    1383 
    1384     buf1p= buf1.rep.elts();
    1385 
    1386     if (lg >= 0 && deggSubLg > 0)
    1387     {
    1388       if (repLengthBuf2 > degfSubLf + 1)
    1389         degfSubLf= repLengthBuf2 - 1;
    1390       tmp= tmin (repLengthBuf1, deggSubLg + 1);
    1391       for (ind= 0; ind < tmp; ind++)
    1392         gp [ind + lg] -= buf1p [ind];
    1393     }
    1394     if (lg < 0)
    1395       break;
    1396 
    1397     buf2p= buf2.rep.elts();
    1398     if (degfSubLf >= 0)
    1399     {
    1400       for (ind= 0; ind < repLengthBuf2; ind++)
    1401         fp [ind + lf] -= buf2p [ind];
    1402     }
    1403   }
    1404 
    1405   return result;
    1406 }
    1407 
    1408 CanonicalForm reverseSubst (const zz_pEX& F, int d, const Variable& alpha)
    1409 {
    1410   Variable y= Variable (2);
    1411   Variable x= Variable (1);
    1412 
    1413   zz_pEX f= F;
    1414   zz_pE *fp= f.rep.elts();
    1415 
    1416   zz_pEX buf;
    1417   zz_pE *bufp;
    1418   CanonicalForm result= 0;
    1419   int i= 0;
    1420   int degf= deg(f);
    1421   int k= 0;
    1422   int degfSubK, repLength, j;
    1423   while (degf >= k)
    1424   {
    1425     degfSubK= degf - k;
    1426     if (degfSubK >= d)
    1427       repLength= d;
    1428     else
    1429       repLength= degfSubK + 1;
    1430 
    1431     buf.rep.SetLength ((long) repLength);
    1432     bufp= buf.rep.elts();
    1433     for (j= 0; j < repLength; j++)
    1434       bufp [j]= fp [j + k];
    1435     buf.normalize();
    1436 
    1437     result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
    1438     i++;
    1439     k= d*i;
    1440   }
    1441 
    1442   return result;
    1443 }
    1444 
    1445 CanonicalForm reverseSubstFp (const zz_pX& F, int d)
    1446 {
    1447   Variable y= Variable (2);
    1448   Variable x= Variable (1);
    1449 
    1450   zz_pX f= F;
    1451   zz_p *fp= f.rep.elts();
    1452 
    1453   zz_pX buf;
    1454   zz_p *bufp;
    1455   CanonicalForm result= 0;
    1456   int i= 0;
    1457   int degf= deg(f);
    1458   int k= 0;
    1459   int degfSubK, repLength, j;
    1460   while (degf >= k)
    1461   {
    1462     degfSubK= degf - k;
    1463     if (degfSubK >= d)
    1464       repLength= d;
    1465     else
    1466       repLength= degfSubK + 1;
    1467 
    1468     buf.rep.SetLength ((long) repLength);
    1469     bufp= buf.rep.elts();
    1470     for (j= 0; j < repLength; j++)
    1471       bufp [j]= fp [j + k];
    1472     buf.normalize();
    1473 
    1474     result += convertNTLzzpX2CF (buf, x)*power (y, i);
    1475     i++;
    1476     k= d*i;
    1477   }
    1478 
    1479   return result;
    1480 }
    1481 
    1482 // assumes input to be reduced mod M and to be an element of Fq not Fp
    1483 CanonicalForm
    1484 mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
    1485                   CanonicalForm& M)
    1486 {
    1487   int d1= degree (F, 1) + degree (G, 1) + 1;
    1488   d1 /= 2;
    1489   d1 += 1;
    1490 
    1491   zz_pX F1, F2;
    1492   kronSubRecipro (F1, F2, F, d1);
    1493   zz_pX G1, G2;
    1494   kronSubRecipro (G1, G2, G, d1);
    1495 
    1496   int k= d1*degree (M);
    1497   MulTrunc (F1, F1, G1, (long) k);
    1498 
    1499   int degtailF= degree (tailcoeff (F), 1);
    1500   int degtailG= degree (tailcoeff (G), 1);
    1501   int taildegF= taildegree (F);
    1502   int taildegG= taildegree (G);
    1503   int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
    1504 
    1505   reverse (F2, F2);
    1506   reverse (G2, G2);
    1507   MulTrunc (F2, F2, G2, b + 1);
    1508   reverse (F2, F2, b);
    1509 
    1510   int d2= tmax (deg (F2)/d1, deg (F1)/d1);
    1511   return reverseSubst (F1, F2, d1, d2);
    1512 }
    1513 
    1514 //Kronecker substitution
    1515 CanonicalForm
    1516 mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
    1517               CanonicalForm& M)
    1518 {
    1519   CanonicalForm A= F;
    1520   CanonicalForm B= G;
    1521 
    1522   int degAx= degree (A, 1);
    1523   int degAy= degree (A, 2);
    1524   int degBx= degree (B, 1);
    1525   int degBy= degree (B, 2);
    1526   int d1= degAx + 1 + degBx;
    1527   int d2= tmax (degAy, degBy);
    1528 
    1529   if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
    1530     return mulMod2NTLFpReci (A, B, M);
    1531 
    1532   zz_pX NTLA= kronSubFp (A, d1);
    1533   zz_pX NTLB= kronSubFp (B, d1);
    1534 
    1535   int k= d1*degree (M);
    1536   MulTrunc (NTLA, NTLA, NTLB, (long) k);
    1537 
    1538   A= reverseSubstFp (NTLA, d1);
    1539 
    1540   return A;
    1541 }
    1542 
    1543 // assumes input to be reduced mod M and to be an element of Fq not Fp
    1544 CanonicalForm
    1545 mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
    1546                   CanonicalForm& M, const Variable& alpha)
    1547 {
    1548   int d1= degree (F, 1) + degree (G, 1) + 1;
    1549   d1 /= 2;
    1550   d1 += 1;
    1551 
    1552   zz_pEX F1, F2;
    1553   kronSubRecipro (F1, F2, F, d1, alpha);
    1554   zz_pEX G1, G2;
    1555   kronSubRecipro (G1, G2, G, d1, alpha);
    1556 
    1557   int k= d1*degree (M);
    1558   MulTrunc (F1, F1, G1, (long) k);
    1559 
    1560   int degtailF= degree (tailcoeff (F), 1);
    1561   int degtailG= degree (tailcoeff (G), 1);
    1562   int taildegF= taildegree (F);
    1563   int taildegG= taildegree (G);
    1564   int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
    1565 
    1566   reverse (F2, F2);
    1567   reverse (G2, G2);
    1568   MulTrunc (F2, F2, G2, b + 1);
    1569   reverse (F2, F2, b);
    1570 
    1571   int d2= tmax (deg (F2)/d1, deg (F1)/d1);
    1572   return reverseSubst (F1, F2, d1, d2, alpha);
    1573 }
    1574 
    1575 #ifdef HAVE_FLINT
    1576 CanonicalForm
    1577 mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
    1578                 CanonicalForm& M);
    1579 #endif
    1580 
    1581 CanonicalForm
    1582 mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
    1583               CanonicalForm& M)
    1584 {
    1585   Variable alpha;
    1586   CanonicalForm A= F;
    1587   CanonicalForm B= G;
    1588 
    1589   if (hasFirstAlgVar (A, alpha) || hasFirstAlgVar (B, alpha))
    1590   {
    1591     int degAx= degree (A, 1);
    1592     int degAy= degree (A, 2);
    1593     int degBx= degree (B, 1);
    1594     int degBy= degree (B, 2);
    1595     int d1= degAx + degBx + 1;
    1596     int d2= tmax (degAy, degBy);
    1597     zz_p::init (getCharacteristic());
    1598     zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
    1599     zz_pE::init (NTLMipo);
    1600 
    1601     int degMipo= degree (getMipo (alpha));
    1602     if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy) &&
    1603         (2*degAy > degree (M)))
    1604       return mulMod2NTLFqReci (A, B, M, alpha);
    1605 
    1606     zz_pEX NTLA= kronSub (A, d1, alpha);
    1607     zz_pEX NTLB= kronSub (B, d1, alpha);
    1608 
    1609     int k= d1*degree (M);
    1610 
    1611     MulTrunc (NTLA, NTLA, NTLB, (long) k);
    1612 
    1613     A= reverseSubst (NTLA, d1, alpha);
    1614 
    1615     return A;
    1616   }
    1617   else
    1618 #ifdef HAVE_FLINT
    1619     return mulMod2FLINTFp (A, B, M);
    1620 #else
    1621     return mulMod2NTLFp (A, B, M);
    1622 #endif
    1623 }
    1624 
    1625 #ifdef HAVE_FLINT
    1626 void
    1627 kronSubRecipro (nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm& A,
    1628                 int d)
    1629 {
    1630   int degAy= degree (A);
    1631   mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
    1632   nmod_poly_init2_preinv (subA1, getCharacteristic(), ninv, d*(degAy + 2));
    1633   nmod_poly_init2_preinv (subA2, getCharacteristic(), ninv, d*(degAy + 2));
    1634 
    1635   nmod_poly_t buf;
    1636 
    1637   int k, kk, j, bufRepLength;
    1638   for (CFIterator i= A; i.hasTerms(); i++)
    1639   {
    1640     convertFacCF2nmod_poly_t (buf, i.coeff());
    1641 
    1642     k= i.exp()*d;
    1643     kk= (degAy - i.exp())*d;
    1644     bufRepLength= (int) nmod_poly_length (buf);
    1645     for (j= 0; j < bufRepLength; j++)
    1646     {
    1647       nmod_poly_set_coeff_ui (subA1, j + k,
    1648                               n_addmod (nmod_poly_get_coeff_ui (subA1, j+k),
    1649                                         nmod_poly_get_coeff_ui (buf, j),
    1650                                         getCharacteristic()
    1651                                        )
    1652                              );
    1653       nmod_poly_set_coeff_ui (subA2, j + kk,
    1654                               n_addmod (nmod_poly_get_coeff_ui (subA2, j + kk),
    1655                                         nmod_poly_get_coeff_ui (buf, j),
    1656                                         getCharacteristic()
    1657                                        )
    1658                              );
    1659     }
    1660     nmod_poly_clear (buf);
    1661   }
    1662   _nmod_poly_normalise (subA1);
    1663   _nmod_poly_normalise (subA2);
    1664 }
    1665 
    1666 void
    1667 kronSubRecipro (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
    1668                 int d)
    1669 {
    1670   int degAy= degree (A);
    1671   fmpz_poly_init2 (subA1, d*(degAy + 2));
    1672   fmpz_poly_init2 (subA2, d*(degAy + 2));
    1673 
    1674   fmpz_poly_t buf;
    1675   fmpz_t coeff1, coeff2;
    1676 
    1677   int k, kk, j, bufRepLength;
    1678   for (CFIterator i= A; i.hasTerms(); i++)
    1679   {
    1680     convertFacCF2Fmpz_poly_t (buf, i.coeff());
    1681 
    1682     k= i.exp()*d;
    1683     kk= (degAy - i.exp())*d;
    1684     bufRepLength= (int) fmpz_poly_length (buf);
    1685     for (j= 0; j < bufRepLength; j++)
    1686     {
    1687       fmpz_poly_get_coeff_fmpz (coeff1, subA1, j+k);
    1688       fmpz_poly_get_coeff_fmpz (coeff2, buf, j);
    1689       fmpz_add (coeff1, coeff1, coeff2);
    1690       fmpz_poly_set_coeff_fmpz (subA1, j + k, coeff1);
    1691       fmpz_poly_get_coeff_fmpz (coeff1, subA2, j + kk);
    1692       fmpz_add (coeff1, coeff1, coeff2);
    1693       fmpz_poly_set_coeff_fmpz (subA2, j + kk, coeff1);
    1694     }
    1695     fmpz_poly_clear (buf);
    1696   }
    1697   fmpz_clear (coeff1);
    1698   fmpz_clear (coeff2);
    1699   _fmpz_poly_normalise (subA1);
    1700   _fmpz_poly_normalise (subA2);
    1701 }
    1702 
    1703 CanonicalForm reverseSubstQ (const fmpz_poly_t F, int d)
    1704 {
    1705   Variable y= Variable (2);
    1706   Variable x= Variable (1);
    1707 
    1708   fmpz_poly_t f;
    1709   fmpz_poly_init (f);
    1710   fmpz_poly_set (f, F);
    1711 
    1712   fmpz_poly_t buf;
    1713   CanonicalForm result= 0;
    1714   int i= 0;
    1715   int degf= fmpz_poly_degree(f);
    1716   int k= 0;
    1717   int degfSubK, repLength, j;
    1718   fmpz_t coeff;
    1719   while (degf >= k)
    1720   {
    1721     degfSubK= degf - k;
    1722     if (degfSubK >= d)
    1723       repLength= d;
    1724     else
    1725       repLength= degfSubK + 1;
    1726 
    1727     fmpz_poly_init2 (buf, repLength);
    1728     fmpz_init (coeff);
    1729     for (j= 0; j < repLength; j++)
    1730     {
    1731       fmpz_poly_get_coeff_fmpz (coeff, f, j + k);
    1732       fmpz_poly_set_coeff_fmpz (buf, j, coeff);
    1733     }
    1734     _fmpz_poly_normalise (buf);
    1735 
    1736     result += convertFmpz_poly_t2FacCF (buf, x)*power (y, i);
    1737     i++;
    1738     k= d*i;
    1739     fmpz_poly_clear (buf);
    1740     fmpz_clear (coeff);
    1741   }
    1742   fmpz_poly_clear (f);
    1743 
    1744   return result;
    1745 }
    1746 
    1747 CanonicalForm
    1748 reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
    1749 {
    1750   Variable y= Variable (2);
    1751   Variable x= Variable (1);
    1752 
    1753   nmod_poly_t f, g;
    1754   mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
    1755   nmod_poly_init_preinv (f, getCharacteristic(), ninv);
    1756   nmod_poly_init_preinv (g, getCharacteristic(), ninv);
    1757   nmod_poly_set (f, F);
    1758   nmod_poly_set (g, G);
    1759   int degf= nmod_poly_degree(f);
    1760   int degg= nmod_poly_degree(g);
    1761 
    1762 
    1763   nmod_poly_t buf1,buf2, buf3;
    1764 
    1765   if (nmod_poly_length (f) < (long) d*(k+1)) //zero padding
    1766     nmod_poly_fit_length (f,(long)d*(k+1));
    1767 
    1768   CanonicalForm result= 0;
    1769   int i= 0;
    1770   int lf= 0;
    1771   int lg= d*k;
    1772   int degfSubLf= degf;
    1773   int deggSubLg= degg-lg;
    1774   int repLengthBuf2, repLengthBuf1, ind, tmp;
    1775   while (degf >= lf || lg >= 0)
    1776   {
    1777     if (degfSubLf >= d)
    1778       repLengthBuf1= d;
    1779     else if (degfSubLf < 0)
    1780       repLengthBuf1= 0;
    1781     else
    1782       repLengthBuf1= degfSubLf + 1;
    1783     nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
    1784 
    1785     for (ind= 0; ind < repLengthBuf1; ind++)
    1786       nmod_poly_set_coeff_ui (buf1, ind, nmod_poly_get_coeff_ui (f, ind+lf));
    1787     _nmod_poly_normalise (buf1);
    1788 
    1789     repLengthBuf1= nmod_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     nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
    1799     for (ind= 0; ind < repLengthBuf2; ind++)
    1800       nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
    1801 
    1802     _nmod_poly_normalise (buf2);
    1803     repLengthBuf2= nmod_poly_length (buf2);
    1804 
    1805     nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
    1806     for (ind= 0; ind < repLengthBuf1; ind++)
    1807       nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
    1808     for (ind= repLengthBuf1; ind < d; ind++)
    1809       nmod_poly_set_coeff_ui (buf3, ind, 0);
    1810     for (ind= 0; ind < repLengthBuf2; ind++)
    1811       nmod_poly_set_coeff_ui (buf3, ind+d, nmod_poly_get_coeff_ui (buf2, ind));
    1812     _nmod_poly_normalise (buf3);
    1813 
    1814     result += convertnmod_poly_t2FacCF (buf3, x)*power (y, i);
    1815     i++;
    1816 
    1817 
    1818     lf= i*d;
    1819     degfSubLf= degf - lf;
    1820 
    1821     lg= d*(k-i);
    1822     deggSubLg= degg - lg;
    1823 
    1824     if (lg >= 0 && deggSubLg > 0)
    1825     {
    1826       if (repLengthBuf2 > degfSubLf + 1)
    1827         degfSubLf= repLengthBuf2 - 1;
    1828       tmp= tmin (repLengthBuf1, deggSubLg + 1);
    1829       for (ind= 0; ind < tmp; ind++)
    1830         nmod_poly_set_coeff_ui (g, ind + lg,
    1831                                 n_submod (nmod_poly_get_coeff_ui (g, ind + lg),
    1832                                           nmod_poly_get_coeff_ui (buf1, ind),
    1833                                           getCharacteristic()
    1834                                          )
    1835                                );
    1836     }
    1837     if (lg < 0)
    1838     {
    1839       nmod_poly_clear (buf1);
    1840       nmod_poly_clear (buf2);
    1841       nmod_poly_clear (buf3);
    1842       break;
    1843     }
    1844     if (degfSubLf >= 0)
    1845     {
    1846       for (ind= 0; ind < repLengthBuf2; ind++)
    1847         nmod_poly_set_coeff_ui (f, ind + lf,
    1848                                 n_submod (nmod_poly_get_coeff_ui (f, ind + lf),
    1849                                           nmod_poly_get_coeff_ui (buf2, ind),
    1850                                           getCharacteristic()
    1851                                          )
    1852                                );
    1853     }
    1854     nmod_poly_clear (buf1);
    1855     nmod_poly_clear (buf2);
    1856     nmod_poly_clear (buf3);
    1857   }
    1858 
    1859   nmod_poly_clear (f);
    1860   nmod_poly_clear (g);
    1861 
    1862   return result;
    1863 }
    1864 
    1865 CanonicalForm
    1866 reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
    1867 {
    1868   Variable y= Variable (2);
    1869   Variable x= Variable (1);
    1870 
    1871   fmpz_poly_t f, g;
    1872   fmpz_poly_init (f);
    1873   fmpz_poly_init (g);
    1874   fmpz_poly_set (f, F);
    1875   fmpz_poly_set (g, G);
    1876   int degf= fmpz_poly_degree(f);
    1877   int degg= fmpz_poly_degree(g);
    1878 
    1879 
    1880   fmpz_poly_t buf1,buf2, buf3;
    1881 
    1882   if (fmpz_poly_length (f) < (long) d*(k+1)) //zero padding
    1883     fmpz_poly_fit_length (f,(long)d*(k+1));
    1884 
    1885   CanonicalForm result= 0;
    1886   int i= 0;
    1887   int lf= 0;
    1888   int lg= d*k;
    1889   int degfSubLf= degf;
    1890   int deggSubLg= degg-lg;
    1891   int repLengthBuf2, repLengthBuf1, ind, tmp;
    1892   fmpz_t tmp1, tmp2;
    1893   while (degf >= lf || lg >= 0)
    1894   {
    1895     if (degfSubLf >= d)
    1896       repLengthBuf1= d;
    1897     else if (degfSubLf < 0)
    1898       repLengthBuf1= 0;
    1899     else
    1900       repLengthBuf1= degfSubLf + 1;
    1901 
    1902     fmpz_poly_init2 (buf1, repLengthBuf1);
    1903 
    1904     for (ind= 0; ind < repLengthBuf1; ind++)
    1905     {
    1906       fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
    1907       fmpz_poly_set_coeff_fmpz (buf1, ind, tmp1);
    1908     }
    1909     _fmpz_poly_normalise (buf1);
    1910 
    1911     repLengthBuf1= fmpz_poly_length (buf1);
    1912 
    1913     if (deggSubLg >= d - 1)
    1914       repLengthBuf2= d - 1;
    1915     else if (deggSubLg < 0)
    1916       repLengthBuf2= 0;
    1917     else
    1918       repLengthBuf2= deggSubLg + 1;
    1919 
    1920     fmpz_poly_init2 (buf2, repLengthBuf2);
    1921 
    1922     for (ind= 0; ind < repLengthBuf2; ind++)
    1923     {
    1924       fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
    1925       fmpz_poly_set_coeff_fmpz (buf2, ind, tmp1);
    1926     }
    1927 
    1928     _fmpz_poly_normalise (buf2);
    1929     repLengthBuf2= fmpz_poly_length (buf2);
    1930 
    1931     fmpz_poly_init2 (buf3, repLengthBuf2 + d);
    1932     for (ind= 0; ind < repLengthBuf1; ind++)
    1933     {
    1934       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))
    1935       fmpz_poly_set_coeff_fmpz (buf3, ind, tmp1);
    1936     }
    1937     for (ind= repLengthBuf1; ind < d; ind++)
    1938       fmpz_poly_set_coeff_ui (buf3, ind, 0);
    1939     for (ind= 0; ind < repLengthBuf2; ind++)
    1940     {
    1941       fmpz_poly_get_coeff_fmpz (tmp1, buf2, ind);
    1942       fmpz_poly_set_coeff_fmpz (buf3, ind + d, tmp1);
    1943     }
    1944     _fmpz_poly_normalise (buf3);
    1945 
    1946     result += convertFmpz_poly_t2FacCF (buf3, x)*power (y, i);
    1947     i++;
    1948 
    1949 
    1950     lf= i*d;
    1951     degfSubLf= degf - lf;
    1952 
    1953     lg= d*(k-i);
    1954     deggSubLg= degg - lg;
    1955 
    1956     if (lg >= 0 && deggSubLg > 0)
    1957     {
    1958       if (repLengthBuf2 > degfSubLf + 1)
    1959         degfSubLf= repLengthBuf2 - 1;
    1960       tmp= tmin (repLengthBuf1, deggSubLg + 1);
    1961       for (ind= 0; ind < tmp; ind++)
    1962       {
    1963         fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
    1964         fmpz_poly_get_coeff_fmpz (tmp2, buf1, ind);
    1965         fmpz_sub (tmp1, tmp1, tmp2);
    1966         fmpz_poly_set_coeff_fmpz (g, ind + lg, tmp1);
    1967       }
    1968     }
    1969     if (lg < 0)
    1970     {
    1971       fmpz_poly_clear (buf1);
    1972       fmpz_poly_clear (buf2);
    1973       fmpz_poly_clear (buf3);
    1974       break;
    1975     }
    1976     if (degfSubLf >= 0)
    1977     {
    1978       for (ind= 0; ind < repLengthBuf2; ind++)
    1979       {
    1980         fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
    1981         fmpz_poly_get_coeff_fmpz (tmp2, buf2, ind);
    1982         fmpz_sub (tmp1, tmp1, tmp2);
    1983         fmpz_poly_set_coeff_fmpz (f, ind + lf, tmp1);
    1984       }
    1985     }
    1986     fmpz_poly_clear (buf1);
    1987     fmpz_poly_clear (buf2);
    1988     fmpz_poly_clear (buf3);
    1989   }
    1990 
    1991   fmpz_poly_clear (f);
    1992   fmpz_poly_clear (g);
    1993   fmpz_clear (tmp1);
    1994   fmpz_clear (tmp2);
    1995 
    1996   return result;
    1997 }
    1998 
    1999 CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
    2000 {
    2001   Variable y= Variable (2);
    2002   Variable x= Variable (1);
    2003 
    2004   nmod_poly_t f;
    2005   mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
    2006   nmod_poly_init_preinv (f, getCharacteristic(), ninv);
    2007   nmod_poly_set (f, F);
    2008 
    2009   nmod_poly_t buf;
    2010   CanonicalForm result= 0;
    2011   int i= 0;
    2012   int degf= nmod_poly_degree(f);
    2013   int k= 0;
    2014   int degfSubK, repLength, j;
    2015   while (degf >= k)
    2016   {
    2017     degfSubK= degf - k;
    2018     if (degfSubK >= d)
    2019       repLength= d;
    2020     else
    2021       repLength= degfSubK + 1;
    2022 
    2023     nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
    2024     for (j= 0; j < repLength; j++)
    2025       nmod_poly_set_coeff_ui (buf, j, nmod_poly_get_coeff_ui (f, j + k));
    2026     _nmod_poly_normalise (buf);
    2027 
    2028     result += convertnmod_poly_t2FacCF (buf, x)*power (y, i);
    2029     i++;
    2030     k= d*i;
    2031     nmod_poly_clear (buf);
    2032   }
    2033   nmod_poly_clear (f);
    2034 
    2035   return result;
    2036 }
    2037 
    2038 CanonicalForm
    2039 mulMod2FLINTFpReci (const CanonicalForm& F, const CanonicalForm& G, const
    2040                     CanonicalForm& M)
    2041 {
    2042   int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
    2043   d1 /= 2;
    2044   d1 += 1;
    2045 
    2046   nmod_poly_t F1, F2;
    2047   mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
    2048   nmod_poly_init_preinv (F1, getCharacteristic(), ninv);
    2049   nmod_poly_init_preinv (F2, getCharacteristic(), ninv);
    2050   kronSubRecipro (F1, F2, F, d1);
    2051 
    2052   nmod_poly_t G1, G2;
    2053   nmod_poly_init_preinv (G1, getCharacteristic(), ninv);
    2054   nmod_poly_init_preinv (G2, getCharacteristic(), ninv);
    2055   kronSubRecipro (G1, G2, G, d1);
    2056 
    2057   int k= d1*degree (M);
    2058   nmod_poly_mullow (F1, F1, G1, (long) k);
    2059 
    2060   int degtailF= degree (tailcoeff (F), 1);;
    2061   int degtailG= degree (tailcoeff (G), 1);
    2062   int taildegF= taildegree (F);
    2063   int taildegG= taildegree (G);
    2064 
    2065   int b= nmod_poly_degree (F2) + nmod_poly_degree (G2) - k - degtailF - degtailG
    2066          + d1*(2+taildegF + taildegG);
    2067   nmod_poly_mulhigh (F2, F2, G2, b);
    2068   nmod_poly_shift_right (F2, F2, b);
    2069   int d2= tmax (nmod_poly_degree (F2)/d1, nmod_poly_degree (F1)/d1);
    2070 
    2071 
    2072   CanonicalForm result= reverseSubst (F1, F2, d1, d2);
    2073 
    2074   nmod_poly_clear (F1);
    2075   nmod_poly_clear (F2);
    2076   nmod_poly_clear (G1);
    2077   nmod_poly_clear (G2);
    2078   return result;
    2079 }
    2080 
    2081 CanonicalForm
    2082 mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
    2083                 CanonicalForm& M)
    2084 {
    2085   CanonicalForm A= F;
    2086   CanonicalForm B= G;
    2087 
    2088   int degAx= degree (A, 1);
    2089   int degAy= degree (A, 2);
    2090   int degBx= degree (B, 1);
    2091   int degBy= degree (B, 2);
    2092   int d1= degAx + 1 + degBx;
    2093   int d2= tmax (degAy, degBy);
    2094 
    2095   if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
    2096     return mulMod2FLINTFpReci (A, B, M);
    2097 
    2098   nmod_poly_t FLINTA, FLINTB;
    2099   mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
    2100   nmod_poly_init_preinv (FLINTA, getCharacteristic(), ninv);
    2101   nmod_poly_init_preinv (FLINTB, getCharacteristic(), ninv);
    2102   kronSubFp (FLINTA, A, d1);
    2103   kronSubFp (FLINTB, B, d1);
    2104 
    2105   int k= d1*degree (M);
    2106   nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
    2107 
    2108   A= reverseSubstFp (FLINTA, d1);
    2109 
    2110   nmod_poly_clear (FLINTA);
    2111   nmod_poly_clear (FLINTB);
    2112   return A;
    2113 }
    2114 
    2115 CanonicalForm
    2116 mulMod2FLINTQReci (const CanonicalForm& F, const CanonicalForm& G, const
    2117                     CanonicalForm& M)
    2118 {
    2119   int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
    2120   d1 /= 2;
    2121   d1 += 1;
    2122 
    2123   fmpz_poly_t F1, F2;
    2124   fmpz_poly_init (F1);
    2125   fmpz_poly_init (F2);
    2126   kronSubRecipro (F1, F2, F, d1);
    2127 
    2128   fmpz_poly_t G1, G2;
    2129   fmpz_poly_init (G1);
    2130   fmpz_poly_init (G2);
    2131   kronSubRecipro (G1, G2, G, d1);
    2132 
    2133   int k= d1*degree (M);
    2134   fmpz_poly_mullow (F1, F1, G1, (long) k);
    2135 
    2136   int degtailF= degree (tailcoeff (F), 1);;
    2137   int degtailG= degree (tailcoeff (G), 1);
    2138   int taildegF= taildegree (F);
    2139   int taildegG= taildegree (G);
    2140 
    2141   int b= fmpz_poly_degree (F2) + fmpz_poly_degree (G2) - k - degtailF - degtailG
    2142          + d1*(2+taildegF + taildegG);
    2143   fmpz_poly_mulhigh_n (F2, F2, G2, b);
    2144   fmpz_poly_shift_right (F2, F2, b);
    2145   int d2= tmax (fmpz_poly_degree (F2)/d1, fmpz_poly_degree (F1)/d1);
    2146 
    2147   CanonicalForm result= reverseSubst (F1, F2, d1, d2);
    2148 
    2149   fmpz_poly_clear (F1);
    2150   fmpz_poly_clear (F2);
    2151   fmpz_poly_clear (G1);
    2152   fmpz_poly_clear (G2);
    2153   return result;
    2154 }
    2155 
    2156 CanonicalForm
    2157 mulMod2FLINTQ (const CanonicalForm& F, const CanonicalForm& G, const
    2158                 CanonicalForm& M)
    2159 {
    2160   CanonicalForm A= F;
    2161   CanonicalForm B= G;
    2162 
    2163   int degAx= degree (A, 1);
    2164   int degBx= degree (B, 1);
    2165   int d1= degAx + 1 + degBx;
    2166 
    2167   CanonicalForm f= bCommonDen (F);
    2168   CanonicalForm g= bCommonDen (G);
    2169   A *= f;
    2170   B *= g;
    2171 
    2172   fmpz_poly_t FLINTA, FLINTB;
    2173   fmpz_poly_init (FLINTA);
    2174   fmpz_poly_init (FLINTB);
    2175   kronSub (FLINTA, A, d1);
    2176   kronSub (FLINTB, B, d1);
    2177   int k= d1*degree (M);
    2178 
    2179   fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
    2180   A= reverseSubstQ (FLINTA, d1);
    2181   fmpz_poly_clear (FLINTA);
    2182   fmpz_poly_clear (FLINTB);
    2183   return A/(f*g);
    2184 }
    2185 
    2186 #endif
    2187 
    2188 CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
    2189                        const CanonicalForm& M)
    2190 {
    2191   if (A.isZero() || B.isZero())
    2192     return 0;
    2193 
    2194   ASSERT (M.isUnivariate(), "M must be univariate");
    2195 
    2196   CanonicalForm F= mod (A, M);
    2197   CanonicalForm G= mod (B, M);
    2198   if (F.inCoeffDomain() || G.inCoeffDomain())
    2199     return F*G;
    2200   Variable y= M.mvar();
    2201   int degF= degree (F, y);
    2202   int degG= degree (G, y);
    2203 
    2204   if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) &&
    2205       (F.level() == G.level()))
    2206   {
    2207     CanonicalForm result= mulNTL (F, G);
    2208     return mod (result, M);
    2209   }
    2210   else if (degF <= 1 && degG <= 1)
    2211   {
    2212     CanonicalForm result= F*G;
    2213     return mod (result, M);
    2214   }
    2215 
    2216   int sizeF= size (F);
    2217   int sizeG= size (G);
    2218 
    2219   int fallBackToNaive= 50;
    2220   if (sizeF < fallBackToNaive || sizeG < fallBackToNaive)
    2221     return mod (F*G, M);
    2222 
    2223 #ifdef HAVE_FLINT
    2224   Variable alpha;
    2225   if (getCharacteristic() == 0 && !hasFirstAlgVar (F, alpha) && ! hasFirstAlgVar (G, alpha))
    2226     return mulMod2FLINTQ (F, G, M);
    2227 #endif
    2228 
    2229   if (getCharacteristic() > 0 && CFFactory::gettype() != GaloisFieldDomain &&
    2230       (((degF-degG) < 50 && degF > degG) || ((degG-degF) < 50 && degF <= degG)))
    2231     return mulMod2NTLFq (F, G, M);
    2232 
    2233   int m= (int) ceil (degree (M)/2.0);
    2234   if (degF >= m || degG >= m)
    2235   {
    2236     CanonicalForm MLo= power (y, m);
    2237     CanonicalForm MHi= power (y, degree (M) - m);
    2238     CanonicalForm F0= mod (F, MLo);
    2239     CanonicalForm F1= div (F, MLo);
    2240     CanonicalForm G0= mod (G, MLo);
    2241     CanonicalForm G1= div (G, MLo);
    2242     CanonicalForm F0G1= mulMod2 (F0, G1, MHi);
    2243     CanonicalForm F1G0= mulMod2 (F1, G0, MHi);
    2244     CanonicalForm F0G0= mulMod2 (F0, G0, M);
    2245     return F0G0 + MLo*(F0G1 + F1G0);
    2246   }
    2247   else
    2248   {
    2249     m= (int) ceil (tmax (degF, degG)/2.0);
    2250     CanonicalForm yToM= power (y, m);
    2251     CanonicalForm F0= mod (F, yToM);
    2252     CanonicalForm F1= div (F, yToM);
    2253     CanonicalForm G0= mod (G, yToM);
    2254     CanonicalForm G1= div (G, yToM);
    2255     CanonicalForm H00= mulMod2 (F0, G0, M);
    2256     CanonicalForm H11= mulMod2 (F1, G1, M);
    2257     CanonicalForm H01= mulMod2 (F0 + F1, G0 + G1, M);
    2258     return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
    2259   }
    2260   DEBOUTLN (cerr, "fatal end in mulMod2");
    2261 }
    2262 
    2263 CanonicalForm mulMod (const CanonicalForm& A, const CanonicalForm& B,
    2264                       const CFList& MOD)
    2265 {
    2266   if (A.isZero() || B.isZero())
    2267     return 0;
    2268 
    2269   if (MOD.length() == 1)
    2270     return mulMod2 (A, B, MOD.getLast());
    2271 
    2272   CanonicalForm M= MOD.getLast();
    2273   CanonicalForm F= mod (A, M);
    2274   CanonicalForm G= mod (B, M);
    2275   if (F.inCoeffDomain() || G.inCoeffDomain())
    2276     return F*G;
    2277   Variable y= M.mvar();
    2278   int degF= degree (F, y);
    2279   int degG= degree (G, y);
    2280 
    2281   if ((degF <= 1 && F.level() <= M.level()) &&
    2282       (degG <= 1 && G.level() <= M.level()))
    2283   {
    2284     CFList buf= MOD;
    2285     buf.removeLast();
    2286     if (degF == 1 && degG == 1)
    2287     {
    2288       CanonicalForm F0= mod (F, y);
    2289       CanonicalForm F1= div (F, y);
    2290       CanonicalForm G0= mod (G, y);
    2291       CanonicalForm G1= div (G, y);
    2292       if (degree (M) > 2)
    2293       {
    2294         CanonicalForm H00= mulMod (F0, G0, buf);
    2295         CanonicalForm H11= mulMod (F1, G1, buf);
    2296         CanonicalForm H01= mulMod (F0 + F1, G0 + G1, buf);
    2297         return H11*y*y + (H01 - H00 - H11)*y + H00;
    2298       }
    2299       else //here degree (M) == 2
    2300       {
    2301         buf.append (y);
    2302         CanonicalForm F0G1= mulMod (F0, G1, buf);
    2303         CanonicalForm F1G0= mulMod (F1, G0, buf);
    2304         CanonicalForm F0G0= mulMod (F0, G0, MOD);
    2305         CanonicalForm result= F0G0 + y*(F0G1 + F1G0);
    2306         return result;
    2307       }
    2308     }
    2309     else if (degF == 1 && degG == 0)
    2310       return mulMod (div (F, y), G, buf)*y + mulMod (mod (F, y), G, buf);
    2311     else if (degF == 0 && degG == 1)
    2312       return mulMod (div (G, y), F, buf)*y + mulMod (mod (G, y), F, buf);
    2313     else
    2314       return mulMod (F, G, buf);
    2315   }
    2316   int m= (int) ceil (degree (M)/2.0);
    2317   if (degF >= m || degG >= m)
    2318   {
    2319     CanonicalForm MLo= power (y, m);
    2320     CanonicalForm MHi= power (y, degree (M) - m);
    2321     CanonicalForm F0= mod (F, MLo);
    2322     CanonicalForm F1= div (F, MLo);
    2323     CanonicalForm G0= mod (G, MLo);
    2324     CanonicalForm G1= div (G, MLo);
    2325     CFList buf= MOD;
    2326     buf.removeLast();
    2327     buf.append (MHi);
    2328     CanonicalForm F0G1= mulMod (F0, G1, buf);
    2329     CanonicalForm F1G0= mulMod (F1, G0, buf);
    2330     CanonicalForm F0G0= mulMod (F0, G0, MOD);
    2331     return F0G0 + MLo*(F0G1 + F1G0);
    2332   }
    2333   else
    2334   {
    2335     m= (int) ceil (tmax (degF, degG)/2.0);
    2336     CanonicalForm yToM= power (y, m);
    2337     CanonicalForm F0= mod (F, yToM);
    2338     CanonicalForm F1= div (F, yToM);
    2339     CanonicalForm G0= mod (G, yToM);
    2340     CanonicalForm G1= div (G, yToM);
    2341     CanonicalForm H00= mulMod (F0, G0, MOD);
    2342     CanonicalForm H11= mulMod (F1, G1, MOD);
    2343     CanonicalForm H01= mulMod (F0 + F1, G0 + G1, MOD);
    2344     return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
    2345   }
    2346   DEBOUTLN (cerr, "fatal end in mulMod");
    2347 }
    2348 
    2349 CanonicalForm prodMod (const CFList& L, const CanonicalForm& M)
    2350 {
    2351   if (L.isEmpty())
    2352     return 1;
    2353   int l= L.length();
    2354   if (l == 1)
    2355     return mod (L.getFirst(), M);
    2356   else if (l == 2) {
    2357     CanonicalForm result= mulMod2 (L.getFirst(), L.getLast(), M);
    2358     return result;
    2359   }
    2360   else
    2361   {
    2362     l /= 2;
    2363     CFList tmp1, tmp2;
    2364     CFListIterator i= L;
    2365     CanonicalForm buf1, buf2;
    2366     for (int j= 1; j <= l; j++, i++)
    2367       tmp1.append (i.getItem());
    2368     tmp2= Difference (L, tmp1);
    2369     buf1= prodMod (tmp1, M);
    2370     buf2= prodMod (tmp2, M);
    2371     CanonicalForm result= mulMod2 (buf1, buf2, M);
    2372     return result;
    2373   }
    2374 }
    2375 
    2376 CanonicalForm prodMod (const CFList& L, const CFList& M)
    2377 {
    2378   if (L.isEmpty())
    2379     return 1;
    2380   else if (L.length() == 1)
    2381     return L.getFirst();
    2382   else if (L.length() == 2)
    2383     return mulMod (L.getFirst(), L.getLast(), M);
    2384   else
    2385   {
    2386     int l= L.length()/2;
    2387     CFListIterator i= L;
    2388     CFList tmp1, tmp2;
    2389     CanonicalForm buf1, buf2;
    2390     for (int j= 1; j <= l; j++, i++)
    2391       tmp1.append (i.getItem());
    2392     tmp2= Difference (L, tmp1);
    2393     buf1= prodMod (tmp1, M);
    2394     buf2= prodMod (tmp2, M);
    2395     return mulMod (buf1, buf2, M);
    2396   }
    2397 }
    2398 
    2399 
    2400 CanonicalForm reverse (const CanonicalForm& F, int d)
    2401 {
    2402   if (d == 0)
    2403     return F;
    2404   CanonicalForm A= F;
    2405   Variable y= Variable (2);
    2406   Variable x= Variable (1);
    2407   if (degree (A, x) > 0)
    2408   {
    2409     A= swapvar (A, x, y);
    2410     CanonicalForm result= 0;
    2411     CFIterator i= A;
    2412     while (d - i.exp() < 0)
    2413       i++;
    2414 
    2415     for (; i.hasTerms() && (d - i.exp() >= 0); i++)
    2416       result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
    2417     return result;
    2418   }
    2419   else
    2420     return A*power (x, d);
    2421 }
    2422 
    2423 CanonicalForm
    2424 newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
    2425 {
    2426   int l= ilog2(n);
    2427 
    2428   CanonicalForm g= mod (F, M)[0] [0];
    2429 
    2430   ASSERT (!g.isZero(), "expected a unit");
    2431 
    2432   Variable alpha;
    2433 
    2434   if (!g.isOne())
    2435     g = 1/g;
    2436   Variable x= Variable (1);
    2437   CanonicalForm result;
    2438   int exp= 0;
    2439   if (n & 1)
    2440   {
    2441     result= g;
    2442     exp= 1;
    2443   }
    2444   CanonicalForm h;
    2445 
    2446   for (int i= 1; i <= l; i++)
    2447   {
    2448     h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
    2449     h= mod (h, power (x, (1 << i)) - 1);
    2450     h= div (h, power (x, (1 << (i - 1))));
    2451     h= mod (h, M);
    2452     g -= power (x, (1 << (i - 1)))*
    2453          mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
    2454 
    2455     if (n & (1 << i))
    2456     {
    2457       if (exp)
    2458       {
    2459         h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
    2460         h= mod (h, power (x, exp + (1 << i)) - 1);
    2461         h= div (h, power (x, exp));
    2462         h= mod (h, M);
    2463         result -= power(x, exp)*mod (mulMod2 (g, h, M),
    2464                                        power (x, (1 << i)));
    2465         exp += (1 << i);
    2466       }
    2467       else
    2468       {
    2469         exp= (1 << i);
    2470         result= g;
    2471       }
    2472     }
    2473   }
    2474 
    2475   return result;
    2476 }
    2477 
    2478 CanonicalForm
    2479 newtonDiv (const CanonicalForm& F, const CanonicalForm& G, const CanonicalForm&
    2480            M)
    2481 {
    2482   ASSERT (getCharacteristic() > 0, "positive characteristic expected");
    2483   ASSERT (CFFactory::gettype() != GaloisFieldDomain, "no GF expected");
    2484 
    2485   CanonicalForm A= mod (F, M);
    2486   CanonicalForm B= mod (G, M);
    2487 
    2488   Variable x= Variable (1);
    2489   int degA= degree (A, x);
    2490   int degB= degree (B, x);
    2491   int m= degA - degB;
    2492   if (m < 0)
    2493     return 0;
    2494 
    2495   Variable v;
    2496   CanonicalForm Q;
    2497   if (degB < 1 || CFFactory::gettype() == GaloisFieldDomain)
    2498   {
    2499     CanonicalForm R;
    2500     divrem2 (A, B, Q, R, M);
    2501   }
    2502   else
    2503   {
    2504     if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
    2505     {
    2506       CanonicalForm R= reverse (A, degA);
    2507       CanonicalForm revB= reverse (B, degB);
    2508       revB= newtonInverse (revB, m + 1, M);
    2509       Q= mulMod2 (R, revB, M);
    2510       Q= mod (Q, power (x, m + 1));
    2511       Q= reverse (Q, m);
    2512     }
    2513     else
    2514     {
    2515       zz_pX mipo= convertFacCF2NTLzzpX (M);
    2516       Variable y= Variable (2);
    2517       zz_pEX NTLA, NTLB;
    2518       NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
    2519       NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
    2520       div (NTLA, NTLA, NTLB);
    2521       Q= convertNTLzz_pEX2CF (NTLA, x, y);
    2522     }
    2523   }
    2524 
    2525   return Q;
    2526 }
    2527 
    2528 void
    2529 newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
    2530               CanonicalForm& R, const CanonicalForm& M)
    2531 {
    2532   CanonicalForm A= mod (F, M);
    2533   CanonicalForm B= mod (G, M);
    2534   Variable x= Variable (1);
    2535   int degA= degree (A, x);
    2536   int degB= degree (B, x);
    2537   int m= degA - degB;
    2538 
    2539   if (m < 0)
    2540   {
    2541     R= A;
    2542     Q= 0;
    2543     return;
    2544   }
    2545 
    2546   Variable v;
    2547   if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
    2548   {
    2549      divrem2 (A, B, Q, R, M);
    2550   }
    2551   else
    2552   {
    2553     if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
    2554     {
    2555       R= reverse (A, degA);
    2556 
    2557       CanonicalForm revB= reverse (B, degB);
    2558       revB= newtonInverse (revB, m + 1, M);
    2559       Q= mulMod2 (R, revB, M);
    2560 
    2561       Q= mod (Q, power (x, m + 1));
    2562       Q= reverse (Q, m);
    2563 
    2564       R= A - mulMod2 (Q, B, M);
    2565     }
    2566     else
    2567     {
    2568       zz_pX mipo= convertFacCF2NTLzzpX (M);
    2569       Variable y= Variable (2);
    2570       zz_pEX NTLA, NTLB;
    2571       NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
    2572       NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
    2573       zz_pEX NTLQ, NTLR;
    2574       DivRem (NTLQ, NTLR, NTLA, NTLB);
    2575       Q= convertNTLzz_pEX2CF (NTLQ, x, y);
    2576       R= convertNTLzz_pEX2CF (NTLR, x, y);
    2577     }
    2578   }
    2579 }
    2580 
    2581 static inline
    2582 CFList split (const CanonicalForm& F, const int m, const Variable& x)
    2583 {
    2584   CanonicalForm A= F;
    2585   CanonicalForm buf= 0;
    2586   bool swap= false;
    2587   if (degree (A, x) <= 0)
    2588     return CFList(A);
    2589   else if (x.level() != A.level())
    2590   {
    2591     swap= true;
    2592     A= swapvar (A, x, A.mvar());
    2593   }
    2594 
    2595   int j= (int) floor ((double) degree (A)/ m);
    2596   CFList result;
    2597   CFIterator i= A;
    2598   for (; j >= 0; j--)
    2599   {
    2600     while (i.hasTerms() && i.exp() - j*m >= 0)
    2601     {
    2602       if (swap)
    2603         buf += i.coeff()*power (A.mvar(), i.exp() - j*m);
    2604       else
    2605         buf += i.coeff()*power (x, i.exp() - j*m);
    2606       i++;
    2607     }
    2608     if (swap)
    2609       result.append (swapvar (buf, x, F.mvar()));
    2610     else
    2611       result.append (buf);
    2612     buf= 0;
    2613   }
    2614   return result;
    2615 }
    2616 
    2617 static inline
    2618 void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
    2619                CanonicalForm& R, const CFList& M);
    2620 
    2621 static inline
    2622 void divrem21 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
    2623                CanonicalForm& R, const CFList& M)
    2624 {
    2625   CanonicalForm A= mod (F, M);
    2626   CanonicalForm B= mod (G, M);
    2627   Variable x= Variable (1);
    2628   int degB= degree (B, x);
    2629   int degA= degree (A, x);
    2630   if (degA < degB)
    2631   {
    2632     Q= 0;
    2633     R= A;
    2634     return;
    2635   }
    2636   ASSERT (2*degB > degA, "expected degree (F, 1) < 2*degree (G, 1)");
    2637   if (degB < 1)
    2638   {
    2639     divrem (A, B, Q, R);
    2640     Q= mod (Q, M);
    2641     R= mod (R, M);
    2642     return;
    2643   }
    2644 
    2645   int m= (int) ceil ((double) (degB + 1)/2.0) + 1;
    2646   CFList splitA= split (A, m, x);
    2647   if (splitA.length() == 3)
    2648     splitA.insert (0);
    2649   if (splitA.length() == 2)
    2650   {
    2651     splitA.insert (0);
    2652     splitA.insert (0);
    2653   }
    2654   if (splitA.length() == 1)
    2655   {
    2656     splitA.insert (0);
    2657     splitA.insert (0);
    2658     splitA.insert (0);
    2659   }
    2660 
    2661   CanonicalForm xToM= power (x, m);
    2662 
    2663   CFListIterator i= splitA;
    2664   CanonicalForm H= i.getItem();
    2665   i++;
    2666   H *= xToM;
    2667   H += i.getItem();
    2668   i++;
    2669   H *= xToM;
    2670   H += i.getItem();
    2671   i++;
    2672 
    2673   divrem32 (H, B, Q, R, M);
    2674 
    2675   CFList splitR= split (R, m, x);
    2676   if (splitR.length() == 1)
    2677     splitR.insert (0);
    2678 
    2679   H= splitR.getFirst();
    2680   H *= xToM;
    2681   H += splitR.getLast();
    2682   H *= xToM;
    2683   H += i.getItem();
    2684 
    2685   CanonicalForm bufQ;
    2686   divrem32 (H, B, bufQ, R, M);
    2687 
    2688   Q *= xToM;
    2689   Q += bufQ;
    2690   return;
    2691 }
    2692 
    2693 static inline
    2694 void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
    2695                CanonicalForm& R, const CFList& M)
    2696 {
    2697   CanonicalForm A= mod (F, M);
    2698   CanonicalForm B= mod (G, M);
    2699   Variable x= Variable (1);
    2700   int degB= degree (B, x);
    2701   int degA= degree (A, x);
    2702   if (degA < degB)
    2703   {
    2704     Q= 0;
    2705     R= A;
    2706     return;
    2707   }
    2708   ASSERT (3*(degB/2) > degA, "expected degree (F, 1) < 3*(degree (G, 1)/2)");
    2709   if (degB < 1)
    2710   {
    2711     divrem (A, B, Q, R);
    2712     Q= mod (Q, M);
    2713     R= mod (R, M);
    2714     return;
    2715   }
    2716   int m= (int) ceil ((double) (degB + 1)/ 2.0);
    2717 
    2718   CFList splitA= split (A, m, x);
    2719   CFList splitB= split (B, m, x);
    2720 
    2721   if (splitA.length() == 2)
    2722   {
    2723     splitA.insert (0);
    2724   }
    2725   if (splitA.length() == 1)
    2726   {
    2727     splitA.insert (0);
    2728     splitA.insert (0);
    2729   }
    2730   CanonicalForm xToM= power (x, m);
    2731 
    2732   CanonicalForm H;
    2733   CFListIterator i= splitA;
    2734   i++;
    2735 
    2736   if (degree (splitA.getFirst(), x) < degree (splitB.getFirst(), x))
    2737   {
    2738     H= splitA.getFirst()*xToM + i.getItem();
    2739     divrem21 (H, splitB.getFirst(), Q, R, M);
    2740   }
    2741   else
    2742   {
    2743     R= splitA.getFirst()*xToM + i.getItem() + splitB.getFirst() -
    2744        splitB.getFirst()*xToM;
    2745     Q= xToM - 1;
    2746   }
    2747 
    2748   H= mulMod (Q, splitB.getLast(), M);
    2749 
    2750   R= R*xToM + splitA.getLast() - H;
    2751 
    2752   while (degree (R, x) >= degB)
    2753   {
    2754     xToM= power (x, degree (R, x) - degB);
    2755     Q += LC (R, x)*xToM;
    2756     R -= mulMod (LC (R, x), B, M)*xToM;
    2757     Q= mod (Q, M);
    2758     R= mod (R, M);
    2759   }
    2760 
    2761   return;
    2762 }
    2763 
    2764 void divrem2 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
    2765               CanonicalForm& R, const CanonicalForm& M)
    2766 {
    2767   CanonicalForm A= mod (F, M);
    2768   CanonicalForm B= mod (G, M);
    2769 
    2770   if (B.inCoeffDomain())
    2771   {
    2772     divrem (A, B, Q, R);
    2773     return;
    2774   }
    2775   if (A.inCoeffDomain() && !B.inCoeffDomain())
    2776   {
    2777     Q= 0;
    2778     R= A;
    2779     return;
    2780   }
    2781 
    2782   if (B.level() < A.level())
    2783   {
    2784     divrem (A, B, Q, R);
    2785     return;
    2786   }
    2787   if (A.level() > B.level())
    2788   {
    2789     R= A;
    2790     Q= 0;
    2791     return;
    2792   }
    2793   if (B.level() == 1 && B.isUnivariate())
    2794   {
    2795     divrem (A, B, Q, R);
    2796     return;
    2797   }
    2798   if (!(B.level() == 1 && B.isUnivariate()) && (A.level() == 1 && A.isUnivariate()))
    2799   {
    2800     Q= 0;
    2801     R= A;
    2802     return;
    2803   }
    2804 
    2805   Variable x= Variable (1);
    2806   int degB= degree (B, x);
    2807   if (degB > degree (A, x))
    2808   {
    2809     Q= 0;
    2810     R= A;
    2811     return;
    2812   }
    2813 
    2814   CFList splitA= split (A, degB, x);
    2815 
    2816   CanonicalForm xToDegB= power (x, degB);
    2817   CanonicalForm H, bufQ;
    2818   Q= 0;
    2819   CFListIterator i= splitA;
    2820   H= i.getItem()*xToDegB;
    2821   i++;
    2822   H += i.getItem();
    2823   CFList buf;
    2824   while (i.hasItem())
    2825   {
    2826     buf= CFList (M);
    2827     divrem21 (H, B, bufQ, R, buf);
    2828     i++;
    2829     if (i.hasItem())
    2830       H= R*xToDegB + i.getItem();
    2831     Q *= xToDegB;
    2832     Q += bufQ;
    2833   }
    2834   return;
    2835 }
    2836 
    2837 void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
    2838              CanonicalForm& R, const CFList& MOD)
    2839 {
    2840   CanonicalForm A= mod (F, MOD);
    2841   CanonicalForm B= mod (G, MOD);
    2842   Variable x= Variable (1);
    2843   int degB= degree (B, x);
    2844   if (degB > degree (A, x))
    2845   {
    2846     Q= 0;
    2847     R= A;
    2848     return;
    2849   }
    2850 
    2851   if (degB <= 0)
    2852   {
    2853     divrem (A, B, Q, R);
    2854     Q= mod (Q, MOD);
    2855     R= mod (R, MOD);
    2856     return;
    2857   }
    2858   CFList splitA= split (A, degB, x);
    2859 
    2860   CanonicalForm xToDegB= power (x, degB);
    2861   CanonicalForm H, bufQ;
    2862   Q= 0;
    2863   CFListIterator i= splitA;
    2864   H= i.getItem()*xToDegB;
    2865   i++;
    2866   H += i.getItem();
    2867   while (i.hasItem())
    2868   {
    2869     divrem21 (H, B, bufQ, R, MOD);
    2870     i++;
    2871     if (i.hasItem())
    2872       H= R*xToDegB + i.getItem();
    2873     Q *= xToDegB;
    2874     Q += bufQ;
    2875   }
    2876   return;
    2877332}
    2878333

  • factory/facHensel.h

    rdeb0f60 r0e2e23  
    44/** @file facHensel.h
    55 *
    6  * This file defines functions for Hensel lifting and modular multiplication.
     6 * This file defines functions for Hensel lifting.
    77 *
    88 * ABSTRACT: function are used for bi- and multivariate factorization over
    … …  
    3030
    3131#ifdef HAVE_NTL
    32 /// multiplication of univariate polys over a finite field using NTL, if we are
    33 /// in GF factory's default multiplication is used.
    34 ///
    35 /// @return @a mulNTL returns F*G
    36 CanonicalForm
    37 mulNTL (const CanonicalForm& F, ///< [in] a univariate poly
    38         const CanonicalForm& G  ///< [in] a univariate poly
    39        );
    40 
    41 /// mod of univariate polys over a finite field using NTL, if we are
    42 /// in GF factory's default mod is used.
    43 ///
    44 /// @return @a modNTL returns F mod G
    45 CanonicalForm
    46 modNTL (const CanonicalForm& F, ///< [in] a univariate poly
    47         const CanonicalForm& G  ///< [in] a univariate poly
    48        );
    49 
    50 /// division of univariate polys over a finite field using NTL, if we are
    51 /// in GF factory's default division is used.
    52 ///
    53 /// @return @a divNTL returns F/G
    54 CanonicalForm
    55 divNTL (const CanonicalForm& F, ///< [in] a univariate poly
    56         const CanonicalForm& G  ///< [in] a univariate poly
    57        );
    58 
    59 /*/// division with remainder of univariate polys over a finite field using NTL,
    60 /// if we are in GF factory's default division with remainder is used.
    61 void
    62 divremNTL (const CanonicalForm& F, ///< [in] a univariate poly
    63            const CanonicalForm& G, ///< [in] a univariate poly
    64            CanonicalForm& Q,       ///< [in,out] dividend
    65            CanonicalForm& R        ///< [in,out] remainder
    66           );*/
    67 
    68 /// division with remainder of @a F by
    69 /// @a G wrt Variable (1) modulo @a M.
    70 ///
    71 /// @return @a Q returns the dividend, @a R returns the remainder.
    72 /// @sa divrem()
    73 void divrem2 (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
    74               const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
    75               CanonicalForm& Q,       ///< [in,out] dividend
    76               CanonicalForm& R,       ///< [in,out] remainder, degree (R, 1) <
    77                                       ///< degree (G, 1)
    78               const CanonicalForm& M  ///< [in] power of Variable (2)
    79              );
    80 
    81 /// division with remainder of @a F by
    82 /// @a G wrt Variable (1) modulo @a MOD.
    83 ///
    84 /// @sa divrem2()
    85 void divrem (
    86            const CanonicalForm& F, ///< [in] multivariate, compressed polynomial
    87            const CanonicalForm& G, ///< [in] multivariate, compressed polynomial
    88            CanonicalForm& Q,       ///< [in,out] dividend
    89            CanonicalForm& R,       ///< [in,out] remainder, degree (R, 1) <
    90                                    ///< degree (G, 1)
    91            const CFList& MOD       ///< [in] only contains powers of
    92                                    ///< Variables of level higher than 1
    93             );
    94 
    95 
    96 /// division with remainder of @a F by
    97 /// @a G wrt Variable (1) modulo @a M using Newton inversion
    98 ///
    99 /// @return @a Q returns the dividend, @a R returns the remainder.
    100 /// @sa divrem2(), newtonDiv()
    101 void
    102 newtonDivrem (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
    103               const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
    104                                       ///< which is monic in Variable (1)
    105               CanonicalForm& Q,       ///< [in,out] dividend
    106               CanonicalForm& R,       ///< [in,out] remainder, degree (R, 1) <
    107                                       ///< degree (G, 1)
    108               const CanonicalForm& M  ///< [in] power of Variable (2)
    109              );
    110 
    111 /// division of @a F by
    112 /// @a G wrt Variable (1) modulo @a M using Newton inversion
    113 ///
    114 /// @return @a newtonDiv returns the dividend
    115 /// @sa divrem2(), newtonDivrem()
    116 CanonicalForm
    117 newtonDiv (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
    118            const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
    119                                    ///< which is monic in Variable (1)
    120            const CanonicalForm& M  ///< [in] power of Variable (2)
    121           );
    122 
    123 /// reduce @a F modulo elements in @a M.
    124 ///
    125 /// @return @a mod returns @a F modulo @a M
    126 CanonicalForm mod (const CanonicalForm& F, ///< [in] compressed polynomial
    127                    const CFList& M         ///< [in] list containing only
    128                                            ///< univariate polynomials
    129                   );
    130 
    131 /// Karatsuba style modular multiplication for bivariate polynomials.
    132 ///
    133 /// @return @a mulMod2 returns @a A * @a B mod @a M.
    134 CanonicalForm
    135 mulMod2 (const CanonicalForm& A, ///< [in] bivariate, compressed polynomial
    136          const CanonicalForm& B, ///< [in] bivariate, compressed polynomial
    137          const CanonicalForm& M  ///< [in] power of Variable (2)
    138         );
    139 
    140 /// Karatsuba style modular multiplication for multivariate polynomials.
    141 ///
    142 /// @return @a mulMod2 returns @a A * @a B mod @a MOD.
    143 CanonicalForm
    144 mulMod (const CanonicalForm& A, ///< [in] multivariate, compressed polynomial
    145         const CanonicalForm& B, ///< [in] multivariate, compressed polynomial
    146         const CFList& MOD       ///< [in] only contains powers of
    147                                 ///< Variables of level higher than 1
    148        );
    149 
    150 /// product of all elements in @a L modulo @a M via divide-and-conquer.
    151 ///
    152 /// @return @a prodMod returns product of all elements in @a L modulo @a M.
    153 CanonicalForm
    154 prodMod (const CFList& L,       ///< [in] contains only bivariate, compressed
    155                                 ///< polynomials
    156          const CanonicalForm& M ///< [in] power of Variable (2)
    157         );
    158 
    159 /// product of all elements in @a L modulo @a M via divide-and-conquer.
    160 ///
    161 /// @return @a prodMod returns product of all elements in @a L modulo @a M.
    162 CanonicalForm
    163 prodMod (const CFList& L, ///< [in] contains multivariate, compressed
    164                           ///< polynomials
    165          const CFList& M  ///< [in] contains only powers of Variables
    166         );
    167 
    16832/// sort a list of polynomials by their degree in @a x.
    16933///
Note: See TracChangeset for help on using the changeset viewer.