[SCM] an open source computer algebra system branch, cleanedupstream, updated. 6125e540ca6d66c307958938a9d53b245507c323
Bernhard R. Link
brlink at debian.org
Tue Apr 24 15:54:47 UTC 2012
The following commit has been merged in the cleanedupstream branch:
commit 5b64ae6261d338d3ec706cdf793dd3def49ea6ce
Author: Martin Lee <martinlee84 at web.de>
Date: Thu Feb 9 16:47:39 2012 +0100
chg: separated multiplication and Hensel lifting functions
diff --git a/factory/GNUmakefile.in b/factory/GNUmakefile.in
index 75ba67e..da2c7c4 100644
--- a/factory/GNUmakefile.in
+++ b/factory/GNUmakefile.in
@@ -180,6 +180,7 @@ basefactorysrc := \
facFqSquarefree.cc \
facHensel.cc \
facIrredTest.cc \
+ facMul.cc \
facSparseHensel.cc \
fieldGCD.cc \
ffops.cc \
@@ -268,6 +269,7 @@ basefactoryincl := \
facFqSquarefree.h \
facHensel.h \
facIrredTest.h \
+ facMul.h \
facSparseHensel.h \
fieldGCD.h \
ffops.h \
diff --git a/factory/facBivar.cc b/factory/facBivar.cc
index eedf9c7..f049c91 100644
--- a/factory/facBivar.cc
+++ b/factory/facBivar.cc
@@ -19,6 +19,7 @@
#include "facFqBivar.h"
#include "facBivar.h"
#include "facHensel.h"
+#include "facMul.h"
#ifdef HAVE_NTL
TIMING_DEFINE_PRINT(fac_uni_factorizer)
diff --git a/factory/facFqBivar.cc b/factory/facFqBivar.cc
index ce3a055..543c475 100644
--- a/factory/facFqBivar.cc
+++ b/factory/facFqBivar.cc
@@ -27,6 +27,7 @@
#include "cf_map_ext.h"
#include "cf_random.h"
#include "facHensel.h"
+#include "facMul.h"
#include "cf_map.h"
#include "cf_gcd_smallp.h"
#include "facFqBivarUtil.h"
diff --git a/factory/facFqBivarUtil.cc b/factory/facFqBivarUtil.cc
index 11ee282..9f48bd4 100644
--- a/factory/facFqBivarUtil.cc
+++ b/factory/facFqBivarUtil.cc
@@ -24,6 +24,7 @@
#include "facFqBivarUtil.h"
#include "cfNewtonPolygon.h"
#include "facHensel.h"
+#include "facMul.h"
void append (CFList& factors1, const CFList& factors2)
diff --git a/factory/facFqFactorize.cc b/factory/facFqFactorize.cc
index 807fa6b..e6351d0 100644
--- a/factory/facFqFactorize.cc
+++ b/factory/facFqFactorize.cc
@@ -29,6 +29,7 @@
#include "cf_map_ext.h"
#include "algext.h"
#include "facSparseHensel.h"
+#include "facMul.h"
#ifdef HAVE_NTL
#include <NTL/ZZ_pEX.h>
diff --git a/factory/facHensel.cc b/factory/facHensel.cc
index 12509b9..bde2c1c 100644
--- a/factory/facHensel.cc
+++ b/factory/facHensel.cc
@@ -3,8 +3,7 @@
\*****************************************************************************/
/** @file facHensel.cc
*
- * This file implements functions to lift factors via Hensel lifting and
- * functions for modular multiplication and division with remainder.
+ * This file implements functions to lift factors via Hensel lifting.
*
* ABSTRACT: Hensel lifting is described in "Efficient Multivariate
* Factorization over Finite Fields" by L. Bernardin & M. Monagon. Division with
@@ -22,6 +21,7 @@
#include "timing.h"
#include "facHensel.h"
+#include "facMul.h"
#include "cf_util.h"
#include "fac_util.h"
#include "cf_algorithm.h"
@@ -31,373 +31,6 @@
#include <NTL/lzz_pEX.h>
#include "NTLconvert.h"
-#ifdef HAVE_FLINT
-#include "FLINTconvert.h"
-#endif
-
-#ifdef HAVE_FLINT
-void kronSub (fmpz_poly_t result, const CanonicalForm& A, int d)
-{
- int degAy= degree (A);
- fmpz_poly_init2 (result, d*(degAy + 1));
- _fmpz_poly_set_length (result, d*(degAy + 1));
- CFIterator j;
- for (CFIterator i= A; i.hasTerms(); i++)
- {
- if (i.coeff().inBaseDomain())
- convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d), i.coeff());
- else
- for (j= i.coeff(); j.hasTerms(); j++)
- convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d+j.exp()),
- j.coeff());
- }
- _fmpz_poly_normalise(result);
-}
-
-
-CanonicalForm
-reverseSubst (const fmpz_poly_t F, int d, const Variable& alpha,
- const CanonicalForm& den)
-{
- Variable x= Variable (1);
-
- CanonicalForm result= 0;
- int i= 0;
- int degf= fmpz_poly_degree (F);
- int k= 0;
- int degfSubK;
- int repLength, j;
- CanonicalForm coeff;
- fmpz* tmp;
- while (degf >= k)
- {
- coeff= 0;
- degfSubK= degf - k;
- if (degfSubK >= d)
- repLength= d;
- else
- repLength= degfSubK + 1;
-
- for (j= 0; j < repLength; j++)
- {
- tmp= fmpz_poly_get_coeff_ptr (F, j+k);
- if (!fmpz_is_zero (tmp))
- {
- CanonicalForm ff= convertFmpz2CF (tmp)/den;
- coeff += ff*power (alpha, j);
- }
- }
- result += coeff*power (x, i);
- i++;
- k= d*i;
- }
- return result;
-}
-
-CanonicalForm
-mulFLINTQa (const CanonicalForm& F, const CanonicalForm& G,
- const Variable& alpha)
-{
- CanonicalForm A= F;
- CanonicalForm B= G;
-
- CanonicalForm denA= bCommonDen (A);
- CanonicalForm denB= bCommonDen (B);
-
- A *= denA;
- B *= denB;
- int degAa= degree (A, alpha);
- int degBa= degree (B, alpha);
- int d= degAa + 1 + degBa;
-
- fmpz_poly_t FLINTA,FLINTB;
- fmpz_poly_init (FLINTA);
- fmpz_poly_init (FLINTB);
- kronSub (FLINTA, A, d);
- kronSub (FLINTB, B, d);
-
- fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
-
- denA *= denB;
- A= reverseSubst (FLINTA, d, alpha, denA);
-
- fmpz_poly_clear (FLINTA);
- fmpz_poly_clear (FLINTB);
- return A;
-}
-
-CanonicalForm
-mulFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
-{
- CanonicalForm A= F;
- CanonicalForm B= G;
-
- CanonicalForm denA= bCommonDen (A);
- CanonicalForm denB= bCommonDen (B);
-
- A *= denA;
- B *= denB;
- fmpz_poly_t FLINTA,FLINTB;
- convertFacCF2Fmpz_poly_t (FLINTA, A);
- convertFacCF2Fmpz_poly_t (FLINTB, B);
- fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
- denA *= denB;
- A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
- A /= denA;
- fmpz_poly_clear (FLINTA);
- fmpz_poly_clear (FLINTB);
-
- return A;
-}
-
-/*CanonicalForm
-mulFLINTQ2 (const CanonicalForm& F, const CanonicalForm& G)
-{
- CanonicalForm A= F;
- CanonicalForm B= G;
-
- fmpq_poly_t FLINTA,FLINTB;
- convertFacCF2Fmpq_poly_t (FLINTA, A);
- convertFacCF2Fmpq_poly_t (FLINTB, B);
-
- fmpq_poly_mul (FLINTA, FLINTA, FLINTB);
- A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
- fmpq_poly_clear (FLINTA);
- fmpq_poly_clear (FLINTB);
- return A;
-}*/
-
-CanonicalForm
-divFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
-{
- CanonicalForm A= F;
- CanonicalForm B= G;
-
- fmpq_poly_t FLINTA,FLINTB;
- convertFacCF2Fmpq_poly_t (FLINTA, A);
- convertFacCF2Fmpq_poly_t (FLINTB, B);
-
- fmpq_poly_div (FLINTA, FLINTA, FLINTB);
- A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
-
- fmpq_poly_clear (FLINTA);
- fmpq_poly_clear (FLINTB);
- return A;
-}
-
-CanonicalForm
-modFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
-{
- CanonicalForm A= F;
- CanonicalForm B= G;
-
- fmpq_poly_t FLINTA,FLINTB;
- convertFacCF2Fmpq_poly_t (FLINTA, A);
- convertFacCF2Fmpq_poly_t (FLINTB, B);
-
- fmpq_poly_rem (FLINTA, FLINTA, FLINTB);
- A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
-
- fmpq_poly_clear (FLINTA);
- fmpq_poly_clear (FLINTB);
- return A;
-}
-
-CanonicalForm
-mulFLINTQaTrunc (const CanonicalForm& F, const CanonicalForm& G,
- const Variable& alpha, int m)
-{
- CanonicalForm A= F;
- CanonicalForm B= G;
-
- CanonicalForm denA= bCommonDen (A);
- CanonicalForm denB= bCommonDen (B);
-
- A *= denA;
- B *= denB;
-
- int degAa= degree (A, alpha);
- int degBa= degree (B, alpha);
- int d= degAa + 1 + degBa;
-
- fmpz_poly_t FLINTA,FLINTB;
- fmpz_poly_init (FLINTA);
- fmpz_poly_init (FLINTB);
- kronSub (FLINTA, A, d);
- kronSub (FLINTB, B, d);
-
- int k= d*m;
- fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, k);
-
- denA *= denB;
- A= reverseSubst (FLINTA, d, alpha, denA);
- fmpz_poly_clear (FLINTA);
- fmpz_poly_clear (FLINTB);
- return A;
-}
-
-CanonicalForm
-mulFLINTQTrunc (const CanonicalForm& F, const CanonicalForm& G, int m)
-{
- if (F.inCoeffDomain() || G.inCoeffDomain())
- return mod (F*G, power (Variable (1), m));
- Variable alpha;
- if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
- return mulFLINTQaTrunc (F, G, alpha, m);
-
- CanonicalForm A= F;
- CanonicalForm B= G;
-
- CanonicalForm denA= bCommonDen (A);
- CanonicalForm denB= bCommonDen (B);
-
- A *= denA;
- B *= denB;
- fmpz_poly_t FLINTA,FLINTB;
- convertFacCF2Fmpz_poly_t (FLINTA, A);
- convertFacCF2Fmpz_poly_t (FLINTB, B);
- fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, m);
- denA *= denB;
- A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
- A /= denA;
- fmpz_poly_clear (FLINTA);
- fmpz_poly_clear (FLINTB);
-
- return A;
-}
-
-CanonicalForm uniReverse (const CanonicalForm& F, int d)
-{
- if (d == 0)
- return F;
- if (F.inCoeffDomain())
- return F*power (Variable (1),d);
- Variable x= Variable (1);
- CanonicalForm result= 0;
- CFIterator i= F;
- while (d - i.exp() < 0)
- i++;
-
- for (; i.hasTerms() && (d - i.exp() >= 0); i++)
- result += i.coeff()*power (x, d - i.exp());
- return result;
-}
-
-CanonicalForm
-newtonInverse (const CanonicalForm& F, const int n)
-{
- int l= ilog2(n);
-
- CanonicalForm g= F [0];
-
- ASSERT (!g.isZero(), "expected a unit");
-
- if (!g.isOne())
- g = 1/g;
- Variable x= Variable (1);
- CanonicalForm result;
- int exp= 0;
- if (n & 1)
- {
- result= g;
- exp= 1;
- }
- CanonicalForm h;
-
- for (int i= 1; i <= l; i++)
- {
- h= mulNTL (g, mod (F, power (x, (1 << i))));
- h= mod (h, power (x, (1 << i)) - 1);
- h= div (h, power (x, (1 << (i - 1))));
- g -= power (x, (1 << (i - 1)))*
- mulFLINTQTrunc (g, h, 1 << (i-1));
-
- if (n & (1 << i))
- {
- if (exp)
- {
- h= mulNTL (result, mod (F, power (x, exp + (1 << i))));
- h= mod (h, power (x, exp + (1 << i)) - 1);
- h= div (h, power (x, exp));
- result -= power(x, exp)*mulFLINTQTrunc (g, h, 1 << i);
- exp += (1 << i);
- }
- else
- {
- exp= (1 << i);
- result= g;
- }
- }
- }
-
- return result;
-}
-
-void
-newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
- CanonicalForm& R)
-{
- CanonicalForm A= F;
- CanonicalForm B= G;
- Variable x= Variable (1);
- int degA= degree (A, x);
- int degB= degree (B, x);
- int m= degA - degB;
-
- if (m < 0)
- {
- R= A;
- Q= 0;
- return;
- }
-
- if (degB <= 1)
- divrem (A, B, Q, R);
- else
- {
- R= uniReverse (A, degA);
-
- CanonicalForm revB= uniReverse (B, degB);
- CanonicalForm buf= revB;
- revB= newtonInverse (revB, m + 1);
- Q= mulFLINTQTrunc (R, revB, m + 1);
- Q= uniReverse (Q, m);
-
- R= A - mulNTL (Q, B);
- }
-}
-
-void
-newtonDiv (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q)
-{
- CanonicalForm A= F;
- CanonicalForm B= G;
- Variable x= Variable (1);
- int degA= degree (A, x);
- int degB= degree (B, x);
- int m= degA - degB;
-
- if (m < 0)
- {
- Q= 0;
- return;
- }
-
- if (degB <= 1)
- Q= div (A, B);
- else
- {
- CanonicalForm R= uniReverse (A, degA);
-
- CanonicalForm revB= uniReverse (B, degB);
- revB= newtonInverse (revB, m + 1);
- Q= mulFLINTQTrunc (R, revB, m + 1);
- Q= uniReverse (Q, m);
- }
-}
-
-#endif
-
static
CFList productsNTL (const CFList& factors, const CanonicalForm& M)
{
@@ -694,2184 +327,6 @@ modularDiophant (const CanonicalForm& f, const CFList& factors,
return result;
}
-CanonicalForm
-mulNTL (const CanonicalForm& F, const CanonicalForm& G)
-{
- if (F.inCoeffDomain() || G.inCoeffDomain() || getCharacteristic() == 0)
- {
- Variable alpha;
-#ifdef HAVE_FLINT
- if ((!F.inCoeffDomain() && !G.inCoeffDomain()) &&
- (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha)))
- {
- CanonicalForm result= mulFLINTQa (F, G, alpha);
- return result;
- }
- else if (!F.inCoeffDomain() && !G.inCoeffDomain())
- return mulFLINTQ (F, G);
-#endif
- return F*G;
- }
- ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
- ASSERT (F.level() == G.level(), "expected polys of same level");
- if (CFFactory::gettype() == GaloisFieldDomain)
- return F*G;
- zz_p::init (getCharacteristic());
- Variable alpha;
- CanonicalForm result;
- if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
- {
- zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
- zz_pE::init (NTLMipo);
- zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
- zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
- mul (NTLF, NTLF, NTLG);
- result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
- }
- else
- {
-#ifdef HAVE_FLINT
- nmod_poly_t FLINTF, FLINTG;
- convertFacCF2nmod_poly_t (FLINTF, F);
- convertFacCF2nmod_poly_t (FLINTG, G);
- nmod_poly_mul (FLINTF, FLINTF, FLINTG);
- result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
- nmod_poly_clear (FLINTF);
- nmod_poly_clear (FLINTG);
-#else
- zz_pX NTLF= convertFacCF2NTLzzpX (F);
- zz_pX NTLG= convertFacCF2NTLzzpX (G);
- mul (NTLF, NTLF, NTLG);
- result= convertNTLzzpX2CF(NTLF, F.mvar());
-#endif
- }
- return result;
-}
-
-CanonicalForm
-modNTL (const CanonicalForm& F, const CanonicalForm& G)
-{
- if (F.inCoeffDomain() && G.isUnivariate())
- return F;
- else if (F.inCoeffDomain() && G.inCoeffDomain())
- return mod (F, G);
- else if (F.isUnivariate() && G.inCoeffDomain())
- return mod (F,G);
-
- if (getCharacteristic() == 0)
- {
-#ifdef HAVE_FLINT
- Variable alpha;
- if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
- return modFLINTQ (F, G);
- else
- {
- CanonicalForm Q, R;
- newtonDivrem (F, G, Q, R);
- return R;
- }
-#else
- return mod (F, G);
-#endif
- }
-
- ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
- ASSERT (F.level() == G.level(), "expected polys of same level");
- if (CFFactory::gettype() == GaloisFieldDomain)
- return mod (F, G);
- zz_p::init (getCharacteristic());
- Variable alpha;
- CanonicalForm result;
- if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
- {
- zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
- zz_pE::init (NTLMipo);
- zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
- zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
- rem (NTLF, NTLF, NTLG);
- result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
- }
- else
- {
-#ifdef HAVE_FLINT
- nmod_poly_t FLINTF, FLINTG;
- convertFacCF2nmod_poly_t (FLINTF, F);
- convertFacCF2nmod_poly_t (FLINTG, G);
- nmod_poly_divrem (FLINTG, FLINTF, FLINTF, FLINTG);
- result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
- nmod_poly_clear (FLINTF);
- nmod_poly_clear (FLINTG);
-#else
- zz_pX NTLF= convertFacCF2NTLzzpX (F);
- zz_pX NTLG= convertFacCF2NTLzzpX (G);
- rem (NTLF, NTLF, NTLG);
- result= convertNTLzzpX2CF(NTLF, F.mvar());
-#endif
- }
- return result;
-}
-
-CanonicalForm
-divNTL (const CanonicalForm& F, const CanonicalForm& G)
-{
- if (F.inCoeffDomain() && G.isUnivariate())
- return F;
- else if (F.inCoeffDomain() && G.inCoeffDomain())
- return div (F, G);
- else if (F.isUnivariate() && G.inCoeffDomain())
- return div (F,G);
-
- if (getCharacteristic() == 0)
- {
-#ifdef HAVE_FLINT
- Variable alpha;
- if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
- return divFLINTQ (F,G);
- else
- {
- CanonicalForm Q;
- newtonDiv (F, G, Q);
- return Q;
- }
-#else
- return div (F, G);
-#endif
- }
-
- ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
- ASSERT (F.level() == G.level(), "expected polys of same level");
- if (CFFactory::gettype() == GaloisFieldDomain)
- return div (F, G);
- zz_p::init (getCharacteristic());
- Variable alpha;
- CanonicalForm result;
- if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
- {
- zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
- zz_pE::init (NTLMipo);
- zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
- zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
- div (NTLF, NTLF, NTLG);
- result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
- }
- else
- {
-#ifdef HAVE_FLINT
- nmod_poly_t FLINTF, FLINTG;
- convertFacCF2nmod_poly_t (FLINTF, F);
- convertFacCF2nmod_poly_t (FLINTG, G);
- nmod_poly_div (FLINTF, FLINTF, FLINTG);
- result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
- nmod_poly_clear (FLINTF);
- nmod_poly_clear (FLINTG);
-#else
- zz_pX NTLF= convertFacCF2NTLzzpX (F);
- zz_pX NTLG= convertFacCF2NTLzzpX (G);
- div (NTLF, NTLF, NTLG);
- result= convertNTLzzpX2CF(NTLF, F.mvar());
-#endif
- }
- return result;
-}
-
-/*
-void
-divremNTL (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
- CanonicalForm& R)
-{
- if (F.inCoeffDomain() && G.isUnivariate())
- {
- R= F;
- Q= 0;
- }
- else if (F.inCoeffDomain() && G.inCoeffDomain())
- {
- divrem (F, G, Q, R);
- return;
- }
- else if (F.isUnivariate() && G.inCoeffDomain())
- {
- divrem (F, G, Q, R);
- return;
- }
-
- ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
- ASSERT (F.level() == G.level(), "expected polys of same level");
- if (CFFactory::gettype() == GaloisFieldDomain)
- {
- divrem (F, G, Q, R);
- return;
- }
- zz_p::init (getCharacteristic());
- Variable alpha;
- CanonicalForm result;
- if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
- {
- zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
- zz_pE::init (NTLMipo);
- zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
- zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
- zz_pEX NTLQ;
- zz_pEX NTLR;
- DivRem (NTLQ, NTLR, NTLF, NTLG);
- Q= convertNTLzz_pEX2CF(NTLQ, F.mvar(), alpha);
- R= convertNTLzz_pEX2CF(NTLR, F.mvar(), alpha);
- return;
- }
- else
- {
- zz_pX NTLF= convertFacCF2NTLzzpX (F);
- zz_pX NTLG= convertFacCF2NTLzzpX (G);
- zz_pX NTLQ;
- zz_pX NTLR;
- DivRem (NTLQ, NTLR, NTLF, NTLG);
- Q= convertNTLzzpX2CF(NTLQ, F.mvar());
- R= convertNTLzzpX2CF(NTLR, F.mvar());
- return;
- }
-}*/
-
-CanonicalForm mod (const CanonicalForm& F, const CFList& M)
-{
- CanonicalForm A= F;
- for (CFListIterator i= M; i.hasItem(); i++)
- A= mod (A, i.getItem());
- return A;
-}
-
-#ifdef HAVE_FLINT
-void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
-{
- int degAy= degree (A);
- nmod_poly_init2 (result, getCharacteristic(), d*(degAy + 1));
-
- nmod_poly_t buf;
-
- int j, k, bufRepLength;
- for (CFIterator i= A; i.hasTerms(); i++)
- {
- convertFacCF2nmod_poly_t (buf, i.coeff());
-
- k= i.exp()*d;
- bufRepLength= (int) nmod_poly_length (buf);
- for (j= 0; j < bufRepLength; j++)
- nmod_poly_set_coeff_ui (result, j + k, nmod_poly_get_coeff_ui (buf, j));
- nmod_poly_clear (buf);
- }
- _nmod_poly_normalise (result);
-}
-
-/*void kronSubQ (fmpz_poly_t result, const CanonicalForm& A, int d)
-{
- int degAy= degree (A);
- fmpz_poly_init2 (result, d*(degAy + 1));
- _fmpz_poly_set_length (result, d*(degAy+1));
-
- CFIterator j;
- for (CFIterator i= A; i.hasTerms(); i++)
- {
- if (i.coeff().inBas
- convertFacCF2Fmpz_poly_t (buf, i.coeff());
-
- int k= i.exp()*d;
- int bufRepLength= (int) fmpz_poly_length (buf);
- for (int j= 0; j < bufRepLength; j++)
- {
- fmpz_poly_get_coeff_fmpz (coeff, buf, j);
- fmpz_poly_set_coeff_fmpz (result, j + k, coeff);
- }
- fmpz_poly_clear (buf);
- }
- fmpz_clear (coeff);
- _fmpz_poly_normalise (result);
-}*/
-
-// A is a bivariate poly over Qa!!!!
-// d2= 2*deg_alpha + 1
-// d1= 2*deg_x*d2+1
-void kronSubQa (fmpq_poly_t result, const CanonicalForm& A, int d1, int d2)
-{
- int degAy= degree (A);
- fmpq_poly_init2 (result, d1*(degAy + 1));
-
- fmpq_poly_t buf;
- fmpq_t coeff;
-
- int k, l, bufRepLength;
- CFIterator j;
- for (CFIterator i= A; i.hasTerms(); i++)
- {
- if (i.coeff().inCoeffDomain())
- {
- k= d1*i.exp();
- convertFacCF2Fmpq_poly_t (buf, i.coeff());
- bufRepLength= (int) fmpq_poly_length(buf);
- for (l= 0; l < bufRepLength; l++)
- {
- fmpq_poly_get_coeff_fmpq (coeff, buf, l);
- fmpq_poly_set_coeff_fmpq (result, l + k, coeff);
- }
- fmpq_poly_clear (buf);
- }
- else
- {
- for (j= i.coeff(); j.hasTerms(); j++)
- {
- k= d1*i.exp();
- k += d2*j.exp();
- convertFacCF2Fmpq_poly_t (buf, j.coeff());
- bufRepLength= (int) fmpq_poly_length(buf);
- for (l= 0; l < bufRepLength; l++)
- {
- fmpq_poly_get_coeff_fmpq (coeff, buf, l);
- fmpq_poly_set_coeff_fmpq (result, k + l, coeff);
- }
- fmpq_poly_clear (buf);
- }
- }
- }
- fmpq_clear (coeff);
- _fmpq_poly_normalise (result);
-}
-#endif
-
-zz_pX kronSubFp (const CanonicalForm& A, int d)
-{
- int degAy= degree (A);
- zz_pX result;
- result.rep.SetLength (d*(degAy + 1));
-
- zz_p *resultp;
- resultp= result.rep.elts();
- zz_pX buf;
- zz_p *bufp;
- int j, k, bufRepLength;
-
- for (CFIterator i= A; i.hasTerms(); i++)
- {
- if (i.coeff().inCoeffDomain())
- buf= convertFacCF2NTLzzpX (i.coeff());
- else
- buf= convertFacCF2NTLzzpX (i.coeff());
-
- k= i.exp()*d;
- bufp= buf.rep.elts();
- bufRepLength= (int) buf.rep.length();
- for (j= 0; j < bufRepLength; j++)
- resultp [j + k]= bufp [j];
- }
- result.normalize();
-
- return result;
-}
-
-zz_pEX kronSub (const CanonicalForm& A, int d, const Variable& alpha)
-{
- int degAy= degree (A);
- zz_pEX result;
- result.rep.SetLength (d*(degAy + 1));
-
- Variable v;
- zz_pE *resultp;
- resultp= result.rep.elts();
- zz_pEX buf1;
- zz_pE *buf1p;
- zz_pX buf2;
- zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
- int j, k, buf1RepLength;
-
- for (CFIterator i= A; i.hasTerms(); i++)
- {
- if (i.coeff().inCoeffDomain())
- {
- buf2= convertFacCF2NTLzzpX (i.coeff());
- buf1= to_zz_pEX (to_zz_pE (buf2));
- }
- else
- buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
-
- k= i.exp()*d;
- buf1p= buf1.rep.elts();
- buf1RepLength= (int) buf1.rep.length();
- for (j= 0; j < buf1RepLength; j++)
- resultp [j + k]= buf1p [j];
- }
- result.normalize();
-
- return result;
-}
-
-void
-kronSubRecipro (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
- const Variable& alpha)
-{
- int degAy= degree (A);
- subA1.rep.SetLength ((long) d*(degAy + 2));
- subA2.rep.SetLength ((long) d*(degAy + 2));
-
- Variable v;
- zz_pE *subA1p;
- zz_pE *subA2p;
- subA1p= subA1.rep.elts();
- subA2p= subA2.rep.elts();
- zz_pEX buf;
- zz_pE *bufp;
- zz_pX buf2;
- zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
- int j, k, kk, bufRepLength;
-
- for (CFIterator i= A; i.hasTerms(); i++)
- {
- if (i.coeff().inCoeffDomain())
- {
- buf2= convertFacCF2NTLzzpX (i.coeff());
- buf= to_zz_pEX (to_zz_pE (buf2));
- }
- else
- buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
-
- k= i.exp()*d;
- kk= (degAy - i.exp())*d;
- bufp= buf.rep.elts();
- bufRepLength= (int) buf.rep.length();
- for (j= 0; j < bufRepLength; j++)
- {
- subA1p [j + k] += bufp [j];
- subA2p [j + kk] += bufp [j];
- }
- }
- subA1.normalize();
- subA2.normalize();
-}
-
-void
-kronSubRecipro (zz_pX& subA1, zz_pX& subA2,const CanonicalForm& A, int d)
-{
- int degAy= degree (A);
- subA1.rep.SetLength ((long) d*(degAy + 2));
- subA2.rep.SetLength ((long) d*(degAy + 2));
-
- zz_p *subA1p;
- zz_p *subA2p;
- subA1p= subA1.rep.elts();
- subA2p= subA2.rep.elts();
- zz_pX buf;
- zz_p *bufp;
- int j, k, kk, bufRepLength;
-
- for (CFIterator i= A; i.hasTerms(); i++)
- {
- buf= convertFacCF2NTLzzpX (i.coeff());
-
- k= i.exp()*d;
- kk= (degAy - i.exp())*d;
- bufp= buf.rep.elts();
- bufRepLength= (int) buf.rep.length();
- for (j= 0; j < bufRepLength; j++)
- {
- subA1p [j + k] += bufp [j];
- subA2p [j + kk] += bufp [j];
- }
- }
- subA1.normalize();
- subA2.normalize();
-}
-
-CanonicalForm
-reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
- const Variable& alpha)
-{
- Variable y= Variable (2);
- Variable x= Variable (1);
-
- zz_pEX f= F;
- zz_pEX g= G;
- int degf= deg(f);
- int degg= deg(g);
-
- zz_pEX buf1;
- zz_pEX buf2;
- zz_pEX buf3;
-
- zz_pE *buf1p;
- zz_pE *buf2p;
- zz_pE *buf3p;
- if (f.rep.length() < (long) d*(k+1)) //zero padding
- f.rep.SetLength ((long)d*(k+1));
-
- zz_pE *gp= g.rep.elts();
- zz_pE *fp= f.rep.elts();
- CanonicalForm result= 0;
- int i= 0;
- int lf= 0;
- int lg= d*k;
- int degfSubLf= degf;
- int deggSubLg= degg-lg;
- int repLengthBuf2, repLengthBuf1, ind, tmp;
- zz_pE zzpEZero= zz_pE();
-
- while (degf >= lf || lg >= 0)
- {
- if (degfSubLf >= d)
- repLengthBuf1= d;
- else if (degfSubLf < 0)
- repLengthBuf1= 0;
- else
- repLengthBuf1= degfSubLf + 1;
- buf1.rep.SetLength((long) repLengthBuf1);
-
- buf1p= buf1.rep.elts();
- for (ind= 0; ind < repLengthBuf1; ind++)
- buf1p [ind]= fp [ind + lf];
- buf1.normalize();
-
- repLengthBuf1= buf1.rep.length();
-
- if (deggSubLg >= d - 1)
- repLengthBuf2= d - 1;
- else if (deggSubLg < 0)
- repLengthBuf2= 0;
- else
- repLengthBuf2= deggSubLg + 1;
-
- buf2.rep.SetLength ((long) repLengthBuf2);
- buf2p= buf2.rep.elts();
- for (ind= 0; ind < repLengthBuf2; ind++)
- buf2p [ind]= gp [ind + lg];
- buf2.normalize();
-
- repLengthBuf2= buf2.rep.length();
-
- buf3.rep.SetLength((long) repLengthBuf2 + d);
- buf3p= buf3.rep.elts();
- buf2p= buf2.rep.elts();
- buf1p= buf1.rep.elts();
- for (ind= 0; ind < repLengthBuf1; ind++)
- buf3p [ind]= buf1p [ind];
- for (ind= repLengthBuf1; ind < d; ind++)
- buf3p [ind]= zzpEZero;
- for (ind= 0; ind < repLengthBuf2; ind++)
- buf3p [ind + d]= buf2p [ind];
- buf3.normalize();
-
- result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
- i++;
-
-
- lf= i*d;
- degfSubLf= degf - lf;
-
- lg= d*(k-i);
- deggSubLg= degg - lg;
-
- buf1p= buf1.rep.elts();
-
- if (lg >= 0 && deggSubLg > 0)
- {
- if (repLengthBuf2 > degfSubLf + 1)
- degfSubLf= repLengthBuf2 - 1;
- tmp= tmin (repLengthBuf1, deggSubLg + 1);
- for (ind= 0; ind < tmp; ind++)
- gp [ind + lg] -= buf1p [ind];
- }
-
- if (lg < 0)
- break;
-
- buf2p= buf2.rep.elts();
- if (degfSubLf >= 0)
- {
- for (ind= 0; ind < repLengthBuf2; ind++)
- fp [ind + lf] -= buf2p [ind];
- }
- }
-
- return result;
-}
-
-CanonicalForm
-reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
-{
- Variable y= Variable (2);
- Variable x= Variable (1);
-
- zz_pX f= F;
- zz_pX g= G;
- int degf= deg(f);
- int degg= deg(g);
-
- zz_pX buf1;
- zz_pX buf2;
- zz_pX buf3;
-
- zz_p *buf1p;
- zz_p *buf2p;
- zz_p *buf3p;
-
- if (f.rep.length() < (long) d*(k+1)) //zero padding
- f.rep.SetLength ((long)d*(k+1));
-
- zz_p *gp= g.rep.elts();
- zz_p *fp= f.rep.elts();
- CanonicalForm result= 0;
- int i= 0;
- int lf= 0;
- int lg= d*k;
- int degfSubLf= degf;
- int deggSubLg= degg-lg;
- int repLengthBuf2, repLengthBuf1, ind, tmp;
- zz_p zzpZero= zz_p();
- while (degf >= lf || lg >= 0)
- {
- if (degfSubLf >= d)
- repLengthBuf1= d;
- else if (degfSubLf < 0)
- repLengthBuf1= 0;
- else
- repLengthBuf1= degfSubLf + 1;
- buf1.rep.SetLength((long) repLengthBuf1);
-
- buf1p= buf1.rep.elts();
- for (ind= 0; ind < repLengthBuf1; ind++)
- buf1p [ind]= fp [ind + lf];
- buf1.normalize();
-
- repLengthBuf1= buf1.rep.length();
-
- if (deggSubLg >= d - 1)
- repLengthBuf2= d - 1;
- else if (deggSubLg < 0)
- repLengthBuf2= 0;
- else
- repLengthBuf2= deggSubLg + 1;
-
- buf2.rep.SetLength ((long) repLengthBuf2);
- buf2p= buf2.rep.elts();
- for (ind= 0; ind < repLengthBuf2; ind++)
- buf2p [ind]= gp [ind + lg];
-
- buf2.normalize();
-
- repLengthBuf2= buf2.rep.length();
-
-
- buf3.rep.SetLength((long) repLengthBuf2 + d);
- buf3p= buf3.rep.elts();
- buf2p= buf2.rep.elts();
- buf1p= buf1.rep.elts();
- for (ind= 0; ind < repLengthBuf1; ind++)
- buf3p [ind]= buf1p [ind];
- for (ind= repLengthBuf1; ind < d; ind++)
- buf3p [ind]= zzpZero;
- for (ind= 0; ind < repLengthBuf2; ind++)
- buf3p [ind + d]= buf2p [ind];
- buf3.normalize();
-
- result += convertNTLzzpX2CF (buf3, x)*power (y, i);
- i++;
-
-
- lf= i*d;
- degfSubLf= degf - lf;
-
- lg= d*(k-i);
- deggSubLg= degg - lg;
-
- buf1p= buf1.rep.elts();
-
- if (lg >= 0 && deggSubLg > 0)
- {
- if (repLengthBuf2 > degfSubLf + 1)
- degfSubLf= repLengthBuf2 - 1;
- tmp= tmin (repLengthBuf1, deggSubLg + 1);
- for (ind= 0; ind < tmp; ind++)
- gp [ind + lg] -= buf1p [ind];
- }
- if (lg < 0)
- break;
-
- buf2p= buf2.rep.elts();
- if (degfSubLf >= 0)
- {
- for (ind= 0; ind < repLengthBuf2; ind++)
- fp [ind + lf] -= buf2p [ind];
- }
- }
-
- return result;
-}
-
-CanonicalForm reverseSubst (const zz_pEX& F, int d, const Variable& alpha)
-{
- Variable y= Variable (2);
- Variable x= Variable (1);
-
- zz_pEX f= F;
- zz_pE *fp= f.rep.elts();
-
- zz_pEX buf;
- zz_pE *bufp;
- CanonicalForm result= 0;
- int i= 0;
- int degf= deg(f);
- int k= 0;
- int degfSubK, repLength, j;
- while (degf >= k)
- {
- degfSubK= degf - k;
- if (degfSubK >= d)
- repLength= d;
- else
- repLength= degfSubK + 1;
-
- buf.rep.SetLength ((long) repLength);
- bufp= buf.rep.elts();
- for (j= 0; j < repLength; j++)
- bufp [j]= fp [j + k];
- buf.normalize();
-
- result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
- i++;
- k= d*i;
- }
-
- return result;
-}
-
-CanonicalForm reverseSubstFp (const zz_pX& F, int d)
-{
- Variable y= Variable (2);
- Variable x= Variable (1);
-
- zz_pX f= F;
- zz_p *fp= f.rep.elts();
-
- zz_pX buf;
- zz_p *bufp;
- CanonicalForm result= 0;
- int i= 0;
- int degf= deg(f);
- int k= 0;
- int degfSubK, repLength, j;
- while (degf >= k)
- {
- degfSubK= degf - k;
- if (degfSubK >= d)
- repLength= d;
- else
- repLength= degfSubK + 1;
-
- buf.rep.SetLength ((long) repLength);
- bufp= buf.rep.elts();
- for (j= 0; j < repLength; j++)
- bufp [j]= fp [j + k];
- buf.normalize();
-
- result += convertNTLzzpX2CF (buf, x)*power (y, i);
- i++;
- k= d*i;
- }
-
- return result;
-}
-
-// assumes input to be reduced mod M and to be an element of Fq not Fp
-CanonicalForm
-mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M)
-{
- int d1= degree (F, 1) + degree (G, 1) + 1;
- d1 /= 2;
- d1 += 1;
-
- zz_pX F1, F2;
- kronSubRecipro (F1, F2, F, d1);
- zz_pX G1, G2;
- kronSubRecipro (G1, G2, G, d1);
-
- int k= d1*degree (M);
- MulTrunc (F1, F1, G1, (long) k);
-
- int degtailF= degree (tailcoeff (F), 1);
- int degtailG= degree (tailcoeff (G), 1);
- int taildegF= taildegree (F);
- int taildegG= taildegree (G);
- int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
-
- reverse (F2, F2);
- reverse (G2, G2);
- MulTrunc (F2, F2, G2, b + 1);
- reverse (F2, F2, b);
-
- int d2= tmax (deg (F2)/d1, deg (F1)/d1);
- return reverseSubst (F1, F2, d1, d2);
-}
-
-//Kronecker substitution
-CanonicalForm
-mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M)
-{
- CanonicalForm A= F;
- CanonicalForm B= G;
-
- int degAx= degree (A, 1);
- int degAy= degree (A, 2);
- int degBx= degree (B, 1);
- int degBy= degree (B, 2);
- int d1= degAx + 1 + degBx;
- int d2= tmax (degAy, degBy);
-
- if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
- return mulMod2NTLFpReci (A, B, M);
-
- zz_pX NTLA= kronSubFp (A, d1);
- zz_pX NTLB= kronSubFp (B, d1);
-
- int k= d1*degree (M);
- MulTrunc (NTLA, NTLA, NTLB, (long) k);
-
- A= reverseSubstFp (NTLA, d1);
-
- return A;
-}
-
-// assumes input to be reduced mod M and to be an element of Fq not Fp
-CanonicalForm
-mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M, const Variable& alpha)
-{
- int d1= degree (F, 1) + degree (G, 1) + 1;
- d1 /= 2;
- d1 += 1;
-
- zz_pEX F1, F2;
- kronSubRecipro (F1, F2, F, d1, alpha);
- zz_pEX G1, G2;
- kronSubRecipro (G1, G2, G, d1, alpha);
-
- int k= d1*degree (M);
- MulTrunc (F1, F1, G1, (long) k);
-
- int degtailF= degree (tailcoeff (F), 1);
- int degtailG= degree (tailcoeff (G), 1);
- int taildegF= taildegree (F);
- int taildegG= taildegree (G);
- int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
-
- reverse (F2, F2);
- reverse (G2, G2);
- MulTrunc (F2, F2, G2, b + 1);
- reverse (F2, F2, b);
-
- int d2= tmax (deg (F2)/d1, deg (F1)/d1);
- return reverseSubst (F1, F2, d1, d2, alpha);
-}
-
-#ifdef HAVE_FLINT
-CanonicalForm
-mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M);
-#endif
-
-CanonicalForm
-mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M)
-{
- Variable alpha;
- CanonicalForm A= F;
- CanonicalForm B= G;
-
- if (hasFirstAlgVar (A, alpha) || hasFirstAlgVar (B, alpha))
- {
- int degAx= degree (A, 1);
- int degAy= degree (A, 2);
- int degBx= degree (B, 1);
- int degBy= degree (B, 2);
- int d1= degAx + degBx + 1;
- int d2= tmax (degAy, degBy);
- zz_p::init (getCharacteristic());
- zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
- zz_pE::init (NTLMipo);
-
- int degMipo= degree (getMipo (alpha));
- if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy) &&
- (2*degAy > degree (M)))
- return mulMod2NTLFqReci (A, B, M, alpha);
-
- zz_pEX NTLA= kronSub (A, d1, alpha);
- zz_pEX NTLB= kronSub (B, d1, alpha);
-
- int k= d1*degree (M);
-
- MulTrunc (NTLA, NTLA, NTLB, (long) k);
-
- A= reverseSubst (NTLA, d1, alpha);
-
- return A;
- }
- else
-#ifdef HAVE_FLINT
- return mulMod2FLINTFp (A, B, M);
-#else
- return mulMod2NTLFp (A, B, M);
-#endif
-}
-
-#ifdef HAVE_FLINT
-void
-kronSubRecipro (nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm& A,
- int d)
-{
- int degAy= degree (A);
- mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
- nmod_poly_init2_preinv (subA1, getCharacteristic(), ninv, d*(degAy + 2));
- nmod_poly_init2_preinv (subA2, getCharacteristic(), ninv, d*(degAy + 2));
-
- nmod_poly_t buf;
-
- int k, kk, j, bufRepLength;
- for (CFIterator i= A; i.hasTerms(); i++)
- {
- convertFacCF2nmod_poly_t (buf, i.coeff());
-
- k= i.exp()*d;
- kk= (degAy - i.exp())*d;
- bufRepLength= (int) nmod_poly_length (buf);
- for (j= 0; j < bufRepLength; j++)
- {
- nmod_poly_set_coeff_ui (subA1, j + k,
- n_addmod (nmod_poly_get_coeff_ui (subA1, j+k),
- nmod_poly_get_coeff_ui (buf, j),
- getCharacteristic()
- )
- );
- nmod_poly_set_coeff_ui (subA2, j + kk,
- n_addmod (nmod_poly_get_coeff_ui (subA2, j + kk),
- nmod_poly_get_coeff_ui (buf, j),
- getCharacteristic()
- )
- );
- }
- nmod_poly_clear (buf);
- }
- _nmod_poly_normalise (subA1);
- _nmod_poly_normalise (subA2);
-}
-
-void
-kronSubRecipro (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
- int d)
-{
- int degAy= degree (A);
- fmpz_poly_init2 (subA1, d*(degAy + 2));
- fmpz_poly_init2 (subA2, d*(degAy + 2));
-
- fmpz_poly_t buf;
- fmpz_t coeff1, coeff2;
-
- int k, kk, j, bufRepLength;
- for (CFIterator i= A; i.hasTerms(); i++)
- {
- convertFacCF2Fmpz_poly_t (buf, i.coeff());
-
- k= i.exp()*d;
- kk= (degAy - i.exp())*d;
- bufRepLength= (int) fmpz_poly_length (buf);
- for (j= 0; j < bufRepLength; j++)
- {
- fmpz_poly_get_coeff_fmpz (coeff1, subA1, j+k);
- fmpz_poly_get_coeff_fmpz (coeff2, buf, j);
- fmpz_add (coeff1, coeff1, coeff2);
- fmpz_poly_set_coeff_fmpz (subA1, j + k, coeff1);
- fmpz_poly_get_coeff_fmpz (coeff1, subA2, j + kk);
- fmpz_add (coeff1, coeff1, coeff2);
- fmpz_poly_set_coeff_fmpz (subA2, j + kk, coeff1);
- }
- fmpz_poly_clear (buf);
- }
- fmpz_clear (coeff1);
- fmpz_clear (coeff2);
- _fmpz_poly_normalise (subA1);
- _fmpz_poly_normalise (subA2);
-}
-
-CanonicalForm reverseSubstQ (const fmpz_poly_t F, int d)
-{
- Variable y= Variable (2);
- Variable x= Variable (1);
-
- fmpz_poly_t f;
- fmpz_poly_init (f);
- fmpz_poly_set (f, F);
-
- fmpz_poly_t buf;
- CanonicalForm result= 0;
- int i= 0;
- int degf= fmpz_poly_degree(f);
- int k= 0;
- int degfSubK, repLength, j;
- fmpz_t coeff;
- while (degf >= k)
- {
- degfSubK= degf - k;
- if (degfSubK >= d)
- repLength= d;
- else
- repLength= degfSubK + 1;
-
- fmpz_poly_init2 (buf, repLength);
- fmpz_init (coeff);
- for (j= 0; j < repLength; j++)
- {
- fmpz_poly_get_coeff_fmpz (coeff, f, j + k);
- fmpz_poly_set_coeff_fmpz (buf, j, coeff);
- }
- _fmpz_poly_normalise (buf);
-
- result += convertFmpz_poly_t2FacCF (buf, x)*power (y, i);
- i++;
- k= d*i;
- fmpz_poly_clear (buf);
- fmpz_clear (coeff);
- }
- fmpz_poly_clear (f);
-
- return result;
-}
-
-CanonicalForm
-reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
-{
- Variable y= Variable (2);
- Variable x= Variable (1);
-
- nmod_poly_t f, g;
- mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
- nmod_poly_init_preinv (f, getCharacteristic(), ninv);
- nmod_poly_init_preinv (g, getCharacteristic(), ninv);
- nmod_poly_set (f, F);
- nmod_poly_set (g, G);
- int degf= nmod_poly_degree(f);
- int degg= nmod_poly_degree(g);
-
-
- nmod_poly_t buf1,buf2, buf3;
-
- if (nmod_poly_length (f) < (long) d*(k+1)) //zero padding
- nmod_poly_fit_length (f,(long)d*(k+1));
-
- CanonicalForm result= 0;
- int i= 0;
- int lf= 0;
- int lg= d*k;
- int degfSubLf= degf;
- int deggSubLg= degg-lg;
- int repLengthBuf2, repLengthBuf1, ind, tmp;
- while (degf >= lf || lg >= 0)
- {
- if (degfSubLf >= d)
- repLengthBuf1= d;
- else if (degfSubLf < 0)
- repLengthBuf1= 0;
- else
- repLengthBuf1= degfSubLf + 1;
- nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
-
- for (ind= 0; ind < repLengthBuf1; ind++)
- nmod_poly_set_coeff_ui (buf1, ind, nmod_poly_get_coeff_ui (f, ind+lf));
- _nmod_poly_normalise (buf1);
-
- repLengthBuf1= nmod_poly_length (buf1);
-
- if (deggSubLg >= d - 1)
- repLengthBuf2= d - 1;
- else if (deggSubLg < 0)
- repLengthBuf2= 0;
- else
- repLengthBuf2= deggSubLg + 1;
-
- nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
- for (ind= 0; ind < repLengthBuf2; ind++)
- nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
-
- _nmod_poly_normalise (buf2);
- repLengthBuf2= nmod_poly_length (buf2);
-
- nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
- for (ind= 0; ind < repLengthBuf1; ind++)
- nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
- for (ind= repLengthBuf1; ind < d; ind++)
- nmod_poly_set_coeff_ui (buf3, ind, 0);
- for (ind= 0; ind < repLengthBuf2; ind++)
- nmod_poly_set_coeff_ui (buf3, ind+d, nmod_poly_get_coeff_ui (buf2, ind));
- _nmod_poly_normalise (buf3);
-
- result += convertnmod_poly_t2FacCF (buf3, x)*power (y, i);
- i++;
-
-
- lf= i*d;
- degfSubLf= degf - lf;
-
- lg= d*(k-i);
- deggSubLg= degg - lg;
-
- if (lg >= 0 && deggSubLg > 0)
- {
- if (repLengthBuf2 > degfSubLf + 1)
- degfSubLf= repLengthBuf2 - 1;
- tmp= tmin (repLengthBuf1, deggSubLg + 1);
- for (ind= 0; ind < tmp; ind++)
- nmod_poly_set_coeff_ui (g, ind + lg,
- n_submod (nmod_poly_get_coeff_ui (g, ind + lg),
- nmod_poly_get_coeff_ui (buf1, ind),
- getCharacteristic()
- )
- );
- }
- if (lg < 0)
- {
- nmod_poly_clear (buf1);
- nmod_poly_clear (buf2);
- nmod_poly_clear (buf3);
- break;
- }
- if (degfSubLf >= 0)
- {
- for (ind= 0; ind < repLengthBuf2; ind++)
- nmod_poly_set_coeff_ui (f, ind + lf,
- n_submod (nmod_poly_get_coeff_ui (f, ind + lf),
- nmod_poly_get_coeff_ui (buf2, ind),
- getCharacteristic()
- )
- );
- }
- nmod_poly_clear (buf1);
- nmod_poly_clear (buf2);
- nmod_poly_clear (buf3);
- }
-
- nmod_poly_clear (f);
- nmod_poly_clear (g);
-
- return result;
-}
-
-CanonicalForm
-reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
-{
- Variable y= Variable (2);
- Variable x= Variable (1);
-
- fmpz_poly_t f, g;
- fmpz_poly_init (f);
- fmpz_poly_init (g);
- fmpz_poly_set (f, F);
- fmpz_poly_set (g, G);
- int degf= fmpz_poly_degree(f);
- int degg= fmpz_poly_degree(g);
-
-
- fmpz_poly_t buf1,buf2, buf3;
-
- if (fmpz_poly_length (f) < (long) d*(k+1)) //zero padding
- fmpz_poly_fit_length (f,(long)d*(k+1));
-
- CanonicalForm result= 0;
- int i= 0;
- int lf= 0;
- int lg= d*k;
- int degfSubLf= degf;
- int deggSubLg= degg-lg;
- int repLengthBuf2, repLengthBuf1, ind, tmp;
- fmpz_t tmp1, tmp2;
- while (degf >= lf || lg >= 0)
- {
- if (degfSubLf >= d)
- repLengthBuf1= d;
- else if (degfSubLf < 0)
- repLengthBuf1= 0;
- else
- repLengthBuf1= degfSubLf + 1;
-
- fmpz_poly_init2 (buf1, repLengthBuf1);
-
- for (ind= 0; ind < repLengthBuf1; ind++)
- {
- fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
- fmpz_poly_set_coeff_fmpz (buf1, ind, tmp1);
- }
- _fmpz_poly_normalise (buf1);
-
- repLengthBuf1= fmpz_poly_length (buf1);
-
- if (deggSubLg >= d - 1)
- repLengthBuf2= d - 1;
- else if (deggSubLg < 0)
- repLengthBuf2= 0;
- else
- repLengthBuf2= deggSubLg + 1;
-
- fmpz_poly_init2 (buf2, repLengthBuf2);
-
- for (ind= 0; ind < repLengthBuf2; ind++)
- {
- fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
- fmpz_poly_set_coeff_fmpz (buf2, ind, tmp1);
- }
-
- _fmpz_poly_normalise (buf2);
- repLengthBuf2= fmpz_poly_length (buf2);
-
- fmpz_poly_init2 (buf3, repLengthBuf2 + d);
- for (ind= 0; ind < repLengthBuf1; ind++)
- {
- 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))
- fmpz_poly_set_coeff_fmpz (buf3, ind, tmp1);
- }
- for (ind= repLengthBuf1; ind < d; ind++)
- fmpz_poly_set_coeff_ui (buf3, ind, 0);
- for (ind= 0; ind < repLengthBuf2; ind++)
- {
- fmpz_poly_get_coeff_fmpz (tmp1, buf2, ind);
- fmpz_poly_set_coeff_fmpz (buf3, ind + d, tmp1);
- }
- _fmpz_poly_normalise (buf3);
-
- result += convertFmpz_poly_t2FacCF (buf3, x)*power (y, i);
- i++;
-
-
- lf= i*d;
- degfSubLf= degf - lf;
-
- lg= d*(k-i);
- deggSubLg= degg - lg;
-
- if (lg >= 0 && deggSubLg > 0)
- {
- if (repLengthBuf2 > degfSubLf + 1)
- degfSubLf= repLengthBuf2 - 1;
- tmp= tmin (repLengthBuf1, deggSubLg + 1);
- for (ind= 0; ind < tmp; ind++)
- {
- fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
- fmpz_poly_get_coeff_fmpz (tmp2, buf1, ind);
- fmpz_sub (tmp1, tmp1, tmp2);
- fmpz_poly_set_coeff_fmpz (g, ind + lg, tmp1);
- }
- }
- if (lg < 0)
- {
- fmpz_poly_clear (buf1);
- fmpz_poly_clear (buf2);
- fmpz_poly_clear (buf3);
- break;
- }
- if (degfSubLf >= 0)
- {
- for (ind= 0; ind < repLengthBuf2; ind++)
- {
- fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
- fmpz_poly_get_coeff_fmpz (tmp2, buf2, ind);
- fmpz_sub (tmp1, tmp1, tmp2);
- fmpz_poly_set_coeff_fmpz (f, ind + lf, tmp1);
- }
- }
- fmpz_poly_clear (buf1);
- fmpz_poly_clear (buf2);
- fmpz_poly_clear (buf3);
- }
-
- fmpz_poly_clear (f);
- fmpz_poly_clear (g);
- fmpz_clear (tmp1);
- fmpz_clear (tmp2);
-
- return result;
-}
-
-CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
-{
- Variable y= Variable (2);
- Variable x= Variable (1);
-
- nmod_poly_t f;
- mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
- nmod_poly_init_preinv (f, getCharacteristic(), ninv);
- nmod_poly_set (f, F);
-
- nmod_poly_t buf;
- CanonicalForm result= 0;
- int i= 0;
- int degf= nmod_poly_degree(f);
- int k= 0;
- int degfSubK, repLength, j;
- while (degf >= k)
- {
- degfSubK= degf - k;
- if (degfSubK >= d)
- repLength= d;
- else
- repLength= degfSubK + 1;
-
- nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
- for (j= 0; j < repLength; j++)
- nmod_poly_set_coeff_ui (buf, j, nmod_poly_get_coeff_ui (f, j + k));
- _nmod_poly_normalise (buf);
-
- result += convertnmod_poly_t2FacCF (buf, x)*power (y, i);
- i++;
- k= d*i;
- nmod_poly_clear (buf);
- }
- nmod_poly_clear (f);
-
- return result;
-}
-
-CanonicalForm
-mulMod2FLINTFpReci (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M)
-{
- int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
- d1 /= 2;
- d1 += 1;
-
- nmod_poly_t F1, F2;
- mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
- nmod_poly_init_preinv (F1, getCharacteristic(), ninv);
- nmod_poly_init_preinv (F2, getCharacteristic(), ninv);
- kronSubRecipro (F1, F2, F, d1);
-
- nmod_poly_t G1, G2;
- nmod_poly_init_preinv (G1, getCharacteristic(), ninv);
- nmod_poly_init_preinv (G2, getCharacteristic(), ninv);
- kronSubRecipro (G1, G2, G, d1);
-
- int k= d1*degree (M);
- nmod_poly_mullow (F1, F1, G1, (long) k);
-
- int degtailF= degree (tailcoeff (F), 1);;
- int degtailG= degree (tailcoeff (G), 1);
- int taildegF= taildegree (F);
- int taildegG= taildegree (G);
-
- int b= nmod_poly_degree (F2) + nmod_poly_degree (G2) - k - degtailF - degtailG
- + d1*(2+taildegF + taildegG);
- nmod_poly_mulhigh (F2, F2, G2, b);
- nmod_poly_shift_right (F2, F2, b);
- int d2= tmax (nmod_poly_degree (F2)/d1, nmod_poly_degree (F1)/d1);
-
-
- CanonicalForm result= reverseSubst (F1, F2, d1, d2);
-
- nmod_poly_clear (F1);
- nmod_poly_clear (F2);
- nmod_poly_clear (G1);
- nmod_poly_clear (G2);
- return result;
-}
-
-CanonicalForm
-mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M)
-{
- CanonicalForm A= F;
- CanonicalForm B= G;
-
- int degAx= degree (A, 1);
- int degAy= degree (A, 2);
- int degBx= degree (B, 1);
- int degBy= degree (B, 2);
- int d1= degAx + 1 + degBx;
- int d2= tmax (degAy, degBy);
-
- if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
- return mulMod2FLINTFpReci (A, B, M);
-
- nmod_poly_t FLINTA, FLINTB;
- mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
- nmod_poly_init_preinv (FLINTA, getCharacteristic(), ninv);
- nmod_poly_init_preinv (FLINTB, getCharacteristic(), ninv);
- kronSubFp (FLINTA, A, d1);
- kronSubFp (FLINTB, B, d1);
-
- int k= d1*degree (M);
- nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
-
- A= reverseSubstFp (FLINTA, d1);
-
- nmod_poly_clear (FLINTA);
- nmod_poly_clear (FLINTB);
- return A;
-}
-
-CanonicalForm
-mulMod2FLINTQReci (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M)
-{
- int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
- d1 /= 2;
- d1 += 1;
-
- fmpz_poly_t F1, F2;
- fmpz_poly_init (F1);
- fmpz_poly_init (F2);
- kronSubRecipro (F1, F2, F, d1);
-
- fmpz_poly_t G1, G2;
- fmpz_poly_init (G1);
- fmpz_poly_init (G2);
- kronSubRecipro (G1, G2, G, d1);
-
- int k= d1*degree (M);
- fmpz_poly_mullow (F1, F1, G1, (long) k);
-
- int degtailF= degree (tailcoeff (F), 1);;
- int degtailG= degree (tailcoeff (G), 1);
- int taildegF= taildegree (F);
- int taildegG= taildegree (G);
-
- int b= fmpz_poly_degree (F2) + fmpz_poly_degree (G2) - k - degtailF - degtailG
- + d1*(2+taildegF + taildegG);
- fmpz_poly_mulhigh_n (F2, F2, G2, b);
- fmpz_poly_shift_right (F2, F2, b);
- int d2= tmax (fmpz_poly_degree (F2)/d1, fmpz_poly_degree (F1)/d1);
-
- CanonicalForm result= reverseSubst (F1, F2, d1, d2);
-
- fmpz_poly_clear (F1);
- fmpz_poly_clear (F2);
- fmpz_poly_clear (G1);
- fmpz_poly_clear (G2);
- return result;
-}
-
-CanonicalForm
-mulMod2FLINTQ (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M)
-{
- CanonicalForm A= F;
- CanonicalForm B= G;
-
- int degAx= degree (A, 1);
- int degBx= degree (B, 1);
- int d1= degAx + 1 + degBx;
-
- CanonicalForm f= bCommonDen (F);
- CanonicalForm g= bCommonDen (G);
- A *= f;
- B *= g;
-
- fmpz_poly_t FLINTA, FLINTB;
- fmpz_poly_init (FLINTA);
- fmpz_poly_init (FLINTB);
- kronSub (FLINTA, A, d1);
- kronSub (FLINTB, B, d1);
- int k= d1*degree (M);
-
- fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
- A= reverseSubstQ (FLINTA, d1);
- fmpz_poly_clear (FLINTA);
- fmpz_poly_clear (FLINTB);
- return A/(f*g);
-}
-
-#endif
-
-CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
- const CanonicalForm& M)
-{
- if (A.isZero() || B.isZero())
- return 0;
-
- ASSERT (M.isUnivariate(), "M must be univariate");
-
- CanonicalForm F= mod (A, M);
- CanonicalForm G= mod (B, M);
- if (F.inCoeffDomain() || G.inCoeffDomain())
- return F*G;
- Variable y= M.mvar();
- int degF= degree (F, y);
- int degG= degree (G, y);
-
- if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) &&
- (F.level() == G.level()))
- {
- CanonicalForm result= mulNTL (F, G);
- return mod (result, M);
- }
- else if (degF <= 1 && degG <= 1)
- {
- CanonicalForm result= F*G;
- return mod (result, M);
- }
-
- int sizeF= size (F);
- int sizeG= size (G);
-
- int fallBackToNaive= 50;
- if (sizeF < fallBackToNaive || sizeG < fallBackToNaive)
- return mod (F*G, M);
-
-#ifdef HAVE_FLINT
- Variable alpha;
- if (getCharacteristic() == 0 && !hasFirstAlgVar (F, alpha) && ! hasFirstAlgVar (G, alpha))
- return mulMod2FLINTQ (F, G, M);
-#endif
-
- if (getCharacteristic() > 0 && CFFactory::gettype() != GaloisFieldDomain &&
- (((degF-degG) < 50 && degF > degG) || ((degG-degF) < 50 && degF <= degG)))
- return mulMod2NTLFq (F, G, M);
-
- int m= (int) ceil (degree (M)/2.0);
- if (degF >= m || degG >= m)
- {
- CanonicalForm MLo= power (y, m);
- CanonicalForm MHi= power (y, degree (M) - m);
- CanonicalForm F0= mod (F, MLo);
- CanonicalForm F1= div (F, MLo);
- CanonicalForm G0= mod (G, MLo);
- CanonicalForm G1= div (G, MLo);
- CanonicalForm F0G1= mulMod2 (F0, G1, MHi);
- CanonicalForm F1G0= mulMod2 (F1, G0, MHi);
- CanonicalForm F0G0= mulMod2 (F0, G0, M);
- return F0G0 + MLo*(F0G1 + F1G0);
- }
- else
- {
- m= (int) ceil (tmax (degF, degG)/2.0);
- CanonicalForm yToM= power (y, m);
- CanonicalForm F0= mod (F, yToM);
- CanonicalForm F1= div (F, yToM);
- CanonicalForm G0= mod (G, yToM);
- CanonicalForm G1= div (G, yToM);
- CanonicalForm H00= mulMod2 (F0, G0, M);
- CanonicalForm H11= mulMod2 (F1, G1, M);
- CanonicalForm H01= mulMod2 (F0 + F1, G0 + G1, M);
- return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
- }
- DEBOUTLN (cerr, "fatal end in mulMod2");
-}
-
-CanonicalForm mulMod (const CanonicalForm& A, const CanonicalForm& B,
- const CFList& MOD)
-{
- if (A.isZero() || B.isZero())
- return 0;
-
- if (MOD.length() == 1)
- return mulMod2 (A, B, MOD.getLast());
-
- CanonicalForm M= MOD.getLast();
- CanonicalForm F= mod (A, M);
- CanonicalForm G= mod (B, M);
- if (F.inCoeffDomain() || G.inCoeffDomain())
- return F*G;
- Variable y= M.mvar();
- int degF= degree (F, y);
- int degG= degree (G, y);
-
- if ((degF <= 1 && F.level() <= M.level()) &&
- (degG <= 1 && G.level() <= M.level()))
- {
- CFList buf= MOD;
- buf.removeLast();
- if (degF == 1 && degG == 1)
- {
- CanonicalForm F0= mod (F, y);
- CanonicalForm F1= div (F, y);
- CanonicalForm G0= mod (G, y);
- CanonicalForm G1= div (G, y);
- if (degree (M) > 2)
- {
- CanonicalForm H00= mulMod (F0, G0, buf);
- CanonicalForm H11= mulMod (F1, G1, buf);
- CanonicalForm H01= mulMod (F0 + F1, G0 + G1, buf);
- return H11*y*y + (H01 - H00 - H11)*y + H00;
- }
- else //here degree (M) == 2
- {
- buf.append (y);
- CanonicalForm F0G1= mulMod (F0, G1, buf);
- CanonicalForm F1G0= mulMod (F1, G0, buf);
- CanonicalForm F0G0= mulMod (F0, G0, MOD);
- CanonicalForm result= F0G0 + y*(F0G1 + F1G0);
- return result;
- }
- }
- else if (degF == 1 && degG == 0)
- return mulMod (div (F, y), G, buf)*y + mulMod (mod (F, y), G, buf);
- else if (degF == 0 && degG == 1)
- return mulMod (div (G, y), F, buf)*y + mulMod (mod (G, y), F, buf);
- else
- return mulMod (F, G, buf);
- }
- int m= (int) ceil (degree (M)/2.0);
- if (degF >= m || degG >= m)
- {
- CanonicalForm MLo= power (y, m);
- CanonicalForm MHi= power (y, degree (M) - m);
- CanonicalForm F0= mod (F, MLo);
- CanonicalForm F1= div (F, MLo);
- CanonicalForm G0= mod (G, MLo);
- CanonicalForm G1= div (G, MLo);
- CFList buf= MOD;
- buf.removeLast();
- buf.append (MHi);
- CanonicalForm F0G1= mulMod (F0, G1, buf);
- CanonicalForm F1G0= mulMod (F1, G0, buf);
- CanonicalForm F0G0= mulMod (F0, G0, MOD);
- return F0G0 + MLo*(F0G1 + F1G0);
- }
- else
- {
- m= (int) ceil (tmax (degF, degG)/2.0);
- CanonicalForm yToM= power (y, m);
- CanonicalForm F0= mod (F, yToM);
- CanonicalForm F1= div (F, yToM);
- CanonicalForm G0= mod (G, yToM);
- CanonicalForm G1= div (G, yToM);
- CanonicalForm H00= mulMod (F0, G0, MOD);
- CanonicalForm H11= mulMod (F1, G1, MOD);
- CanonicalForm H01= mulMod (F0 + F1, G0 + G1, MOD);
- return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
- }
- DEBOUTLN (cerr, "fatal end in mulMod");
-}
-
-CanonicalForm prodMod (const CFList& L, const CanonicalForm& M)
-{
- if (L.isEmpty())
- return 1;
- int l= L.length();
- if (l == 1)
- return mod (L.getFirst(), M);
- else if (l == 2) {
- CanonicalForm result= mulMod2 (L.getFirst(), L.getLast(), M);
- return result;
- }
- else
- {
- l /= 2;
- CFList tmp1, tmp2;
- CFListIterator i= L;
- CanonicalForm buf1, buf2;
- for (int j= 1; j <= l; j++, i++)
- tmp1.append (i.getItem());
- tmp2= Difference (L, tmp1);
- buf1= prodMod (tmp1, M);
- buf2= prodMod (tmp2, M);
- CanonicalForm result= mulMod2 (buf1, buf2, M);
- return result;
- }
-}
-
-CanonicalForm prodMod (const CFList& L, const CFList& M)
-{
- if (L.isEmpty())
- return 1;
- else if (L.length() == 1)
- return L.getFirst();
- else if (L.length() == 2)
- return mulMod (L.getFirst(), L.getLast(), M);
- else
- {
- int l= L.length()/2;
- CFListIterator i= L;
- CFList tmp1, tmp2;
- CanonicalForm buf1, buf2;
- for (int j= 1; j <= l; j++, i++)
- tmp1.append (i.getItem());
- tmp2= Difference (L, tmp1);
- buf1= prodMod (tmp1, M);
- buf2= prodMod (tmp2, M);
- return mulMod (buf1, buf2, M);
- }
-}
-
-
-CanonicalForm reverse (const CanonicalForm& F, int d)
-{
- if (d == 0)
- return F;
- CanonicalForm A= F;
- Variable y= Variable (2);
- Variable x= Variable (1);
- if (degree (A, x) > 0)
- {
- A= swapvar (A, x, y);
- CanonicalForm result= 0;
- CFIterator i= A;
- while (d - i.exp() < 0)
- i++;
-
- for (; i.hasTerms() && (d - i.exp() >= 0); i++)
- result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
- return result;
- }
- else
- return A*power (x, d);
-}
-
-CanonicalForm
-newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
-{
- int l= ilog2(n);
-
- CanonicalForm g= mod (F, M)[0] [0];
-
- ASSERT (!g.isZero(), "expected a unit");
-
- Variable alpha;
-
- if (!g.isOne())
- g = 1/g;
- Variable x= Variable (1);
- CanonicalForm result;
- int exp= 0;
- if (n & 1)
- {
- result= g;
- exp= 1;
- }
- CanonicalForm h;
-
- for (int i= 1; i <= l; i++)
- {
- h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
- h= mod (h, power (x, (1 << i)) - 1);
- h= div (h, power (x, (1 << (i - 1))));
- h= mod (h, M);
- g -= power (x, (1 << (i - 1)))*
- mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
-
- if (n & (1 << i))
- {
- if (exp)
- {
- h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
- h= mod (h, power (x, exp + (1 << i)) - 1);
- h= div (h, power (x, exp));
- h= mod (h, M);
- result -= power(x, exp)*mod (mulMod2 (g, h, M),
- power (x, (1 << i)));
- exp += (1 << i);
- }
- else
- {
- exp= (1 << i);
- result= g;
- }
- }
- }
-
- return result;
-}
-
-CanonicalForm
-newtonDiv (const CanonicalForm& F, const CanonicalForm& G, const CanonicalForm&
- M)
-{
- ASSERT (getCharacteristic() > 0, "positive characteristic expected");
- ASSERT (CFFactory::gettype() != GaloisFieldDomain, "no GF expected");
-
- CanonicalForm A= mod (F, M);
- CanonicalForm B= mod (G, M);
-
- Variable x= Variable (1);
- int degA= degree (A, x);
- int degB= degree (B, x);
- int m= degA - degB;
- if (m < 0)
- return 0;
-
- Variable v;
- CanonicalForm Q;
- if (degB < 1 || CFFactory::gettype() == GaloisFieldDomain)
- {
- CanonicalForm R;
- divrem2 (A, B, Q, R, M);
- }
- else
- {
- if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
- {
- CanonicalForm R= reverse (A, degA);
- CanonicalForm revB= reverse (B, degB);
- revB= newtonInverse (revB, m + 1, M);
- Q= mulMod2 (R, revB, M);
- Q= mod (Q, power (x, m + 1));
- Q= reverse (Q, m);
- }
- else
- {
- zz_pX mipo= convertFacCF2NTLzzpX (M);
- Variable y= Variable (2);
- zz_pEX NTLA, NTLB;
- NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
- NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
- div (NTLA, NTLA, NTLB);
- Q= convertNTLzz_pEX2CF (NTLA, x, y);
- }
- }
-
- return Q;
-}
-
-void
-newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
- CanonicalForm& R, const CanonicalForm& M)
-{
- CanonicalForm A= mod (F, M);
- CanonicalForm B= mod (G, M);
- Variable x= Variable (1);
- int degA= degree (A, x);
- int degB= degree (B, x);
- int m= degA - degB;
-
- if (m < 0)
- {
- R= A;
- Q= 0;
- return;
- }
-
- Variable v;
- if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
- {
- divrem2 (A, B, Q, R, M);
- }
- else
- {
- if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
- {
- R= reverse (A, degA);
-
- CanonicalForm revB= reverse (B, degB);
- revB= newtonInverse (revB, m + 1, M);
- Q= mulMod2 (R, revB, M);
-
- Q= mod (Q, power (x, m + 1));
- Q= reverse (Q, m);
-
- R= A - mulMod2 (Q, B, M);
- }
- else
- {
- zz_pX mipo= convertFacCF2NTLzzpX (M);
- Variable y= Variable (2);
- zz_pEX NTLA, NTLB;
- NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
- NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
- zz_pEX NTLQ, NTLR;
- DivRem (NTLQ, NTLR, NTLA, NTLB);
- Q= convertNTLzz_pEX2CF (NTLQ, x, y);
- R= convertNTLzz_pEX2CF (NTLR, x, y);
- }
- }
-}
-
-static inline
-CFList split (const CanonicalForm& F, const int m, const Variable& x)
-{
- CanonicalForm A= F;
- CanonicalForm buf= 0;
- bool swap= false;
- if (degree (A, x) <= 0)
- return CFList(A);
- else if (x.level() != A.level())
- {
- swap= true;
- A= swapvar (A, x, A.mvar());
- }
-
- int j= (int) floor ((double) degree (A)/ m);
- CFList result;
- CFIterator i= A;
- for (; j >= 0; j--)
- {
- while (i.hasTerms() && i.exp() - j*m >= 0)
- {
- if (swap)
- buf += i.coeff()*power (A.mvar(), i.exp() - j*m);
- else
- buf += i.coeff()*power (x, i.exp() - j*m);
- i++;
- }
- if (swap)
- result.append (swapvar (buf, x, F.mvar()));
- else
- result.append (buf);
- buf= 0;
- }
- return result;
-}
-
-static inline
-void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
- CanonicalForm& R, const CFList& M);
-
-static inline
-void divrem21 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
- CanonicalForm& R, const CFList& M)
-{
- CanonicalForm A= mod (F, M);
- CanonicalForm B= mod (G, M);
- Variable x= Variable (1);
- int degB= degree (B, x);
- int degA= degree (A, x);
- if (degA < degB)
- {
- Q= 0;
- R= A;
- return;
- }
- ASSERT (2*degB > degA, "expected degree (F, 1) < 2*degree (G, 1)");
- if (degB < 1)
- {
- divrem (A, B, Q, R);
- Q= mod (Q, M);
- R= mod (R, M);
- return;
- }
-
- int m= (int) ceil ((double) (degB + 1)/2.0) + 1;
- CFList splitA= split (A, m, x);
- if (splitA.length() == 3)
- splitA.insert (0);
- if (splitA.length() == 2)
- {
- splitA.insert (0);
- splitA.insert (0);
- }
- if (splitA.length() == 1)
- {
- splitA.insert (0);
- splitA.insert (0);
- splitA.insert (0);
- }
-
- CanonicalForm xToM= power (x, m);
-
- CFListIterator i= splitA;
- CanonicalForm H= i.getItem();
- i++;
- H *= xToM;
- H += i.getItem();
- i++;
- H *= xToM;
- H += i.getItem();
- i++;
-
- divrem32 (H, B, Q, R, M);
-
- CFList splitR= split (R, m, x);
- if (splitR.length() == 1)
- splitR.insert (0);
-
- H= splitR.getFirst();
- H *= xToM;
- H += splitR.getLast();
- H *= xToM;
- H += i.getItem();
-
- CanonicalForm bufQ;
- divrem32 (H, B, bufQ, R, M);
-
- Q *= xToM;
- Q += bufQ;
- return;
-}
-
-static inline
-void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
- CanonicalForm& R, const CFList& M)
-{
- CanonicalForm A= mod (F, M);
- CanonicalForm B= mod (G, M);
- Variable x= Variable (1);
- int degB= degree (B, x);
- int degA= degree (A, x);
- if (degA < degB)
- {
- Q= 0;
- R= A;
- return;
- }
- ASSERT (3*(degB/2) > degA, "expected degree (F, 1) < 3*(degree (G, 1)/2)");
- if (degB < 1)
- {
- divrem (A, B, Q, R);
- Q= mod (Q, M);
- R= mod (R, M);
- return;
- }
- int m= (int) ceil ((double) (degB + 1)/ 2.0);
-
- CFList splitA= split (A, m, x);
- CFList splitB= split (B, m, x);
-
- if (splitA.length() == 2)
- {
- splitA.insert (0);
- }
- if (splitA.length() == 1)
- {
- splitA.insert (0);
- splitA.insert (0);
- }
- CanonicalForm xToM= power (x, m);
-
- CanonicalForm H;
- CFListIterator i= splitA;
- i++;
-
- if (degree (splitA.getFirst(), x) < degree (splitB.getFirst(), x))
- {
- H= splitA.getFirst()*xToM + i.getItem();
- divrem21 (H, splitB.getFirst(), Q, R, M);
- }
- else
- {
- R= splitA.getFirst()*xToM + i.getItem() + splitB.getFirst() -
- splitB.getFirst()*xToM;
- Q= xToM - 1;
- }
-
- H= mulMod (Q, splitB.getLast(), M);
-
- R= R*xToM + splitA.getLast() - H;
-
- while (degree (R, x) >= degB)
- {
- xToM= power (x, degree (R, x) - degB);
- Q += LC (R, x)*xToM;
- R -= mulMod (LC (R, x), B, M)*xToM;
- Q= mod (Q, M);
- R= mod (R, M);
- }
-
- return;
-}
-
-void divrem2 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
- CanonicalForm& R, const CanonicalForm& M)
-{
- CanonicalForm A= mod (F, M);
- CanonicalForm B= mod (G, M);
-
- if (B.inCoeffDomain())
- {
- divrem (A, B, Q, R);
- return;
- }
- if (A.inCoeffDomain() && !B.inCoeffDomain())
- {
- Q= 0;
- R= A;
- return;
- }
-
- if (B.level() < A.level())
- {
- divrem (A, B, Q, R);
- return;
- }
- if (A.level() > B.level())
- {
- R= A;
- Q= 0;
- return;
- }
- if (B.level() == 1 && B.isUnivariate())
- {
- divrem (A, B, Q, R);
- return;
- }
- if (!(B.level() == 1 && B.isUnivariate()) && (A.level() == 1 && A.isUnivariate()))
- {
- Q= 0;
- R= A;
- return;
- }
-
- Variable x= Variable (1);
- int degB= degree (B, x);
- if (degB > degree (A, x))
- {
- Q= 0;
- R= A;
- return;
- }
-
- CFList splitA= split (A, degB, x);
-
- CanonicalForm xToDegB= power (x, degB);
- CanonicalForm H, bufQ;
- Q= 0;
- CFListIterator i= splitA;
- H= i.getItem()*xToDegB;
- i++;
- H += i.getItem();
- CFList buf;
- while (i.hasItem())
- {
- buf= CFList (M);
- divrem21 (H, B, bufQ, R, buf);
- i++;
- if (i.hasItem())
- H= R*xToDegB + i.getItem();
- Q *= xToDegB;
- Q += bufQ;
- }
- return;
-}
-
-void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
- CanonicalForm& R, const CFList& MOD)
-{
- CanonicalForm A= mod (F, MOD);
- CanonicalForm B= mod (G, MOD);
- Variable x= Variable (1);
- int degB= degree (B, x);
- if (degB > degree (A, x))
- {
- Q= 0;
- R= A;
- return;
- }
-
- if (degB <= 0)
- {
- divrem (A, B, Q, R);
- Q= mod (Q, MOD);
- R= mod (R, MOD);
- return;
- }
- CFList splitA= split (A, degB, x);
-
- CanonicalForm xToDegB= power (x, degB);
- CanonicalForm H, bufQ;
- Q= 0;
- CFListIterator i= splitA;
- H= i.getItem()*xToDegB;
- i++;
- H += i.getItem();
- while (i.hasItem())
- {
- divrem21 (H, B, bufQ, R, MOD);
- i++;
- if (i.hasItem())
- H= R*xToDegB + i.getItem();
- Q *= xToDegB;
- Q += bufQ;
- }
- return;
-}
-
void sortList (CFList& list, const Variable& x)
{
int l= 1;
diff --git a/factory/facHensel.h b/factory/facHensel.h
index bf495d4..e3bf585 100644
--- a/factory/facHensel.h
+++ b/factory/facHensel.h
@@ -3,7 +3,7 @@
\*****************************************************************************/
/** @file facHensel.h
*
- * This file defines functions for Hensel lifting and modular multiplication.
+ * This file defines functions for Hensel lifting.
*
* ABSTRACT: function are used for bi- and multivariate factorization over
* finite fields
@@ -26,142 +26,6 @@
#include "templates/ftmpl_functions.h"
#include "algext.h"
-/// multiplication of univariate polys over a finite field using NTL, if we are
-/// in GF factory's default multiplication is used.
-///
-/// @return @a mulNTL returns F*G
-CanonicalForm
-mulNTL (const CanonicalForm& F, ///< [in] a univariate poly
- const CanonicalForm& G ///< [in] a univariate poly
- );
-
-/// mod of univariate polys over a finite field using NTL, if we are
-/// in GF factory's default mod is used.
-///
-/// @return @a modNTL returns F mod G
-CanonicalForm
-modNTL (const CanonicalForm& F, ///< [in] a univariate poly
- const CanonicalForm& G ///< [in] a univariate poly
- );
-
-/// division of univariate polys over a finite field using NTL, if we are
-/// in GF factory's default division is used.
-///
-/// @return @a divNTL returns F/G
-CanonicalForm
-divNTL (const CanonicalForm& F, ///< [in] a univariate poly
- const CanonicalForm& G ///< [in] a univariate poly
- );
-
-/*/// division with remainder of univariate polys over a finite field using NTL,
-/// if we are in GF factory's default division with remainder is used.
-void
-divremNTL (const CanonicalForm& F, ///< [in] a univariate poly
- const CanonicalForm& G, ///< [in] a univariate poly
- CanonicalForm& Q, ///< [in,out] dividend
- CanonicalForm& R ///< [in,out] remainder
- );*/
-
-/// division with remainder of @a F by
-/// @a G wrt Variable (1) modulo @a M.
-///
-/// @return @a Q returns the dividend, @a R returns the remainder.
-/// @sa divrem()
-void divrem2 (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
- const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
- CanonicalForm& Q, ///< [in,out] dividend
- CanonicalForm& R, ///< [in,out] remainder, degree (R, 1) <
- ///< degree (G, 1)
- const CanonicalForm& M ///< [in] power of Variable (2)
- );
-
-/// division with remainder of @a F by
-/// @a G wrt Variable (1) modulo @a MOD.
-///
-/// @sa divrem2()
-void divrem (
- const CanonicalForm& F, ///< [in] multivariate, compressed polynomial
- const CanonicalForm& G, ///< [in] multivariate, compressed polynomial
- CanonicalForm& Q, ///< [in,out] dividend
- CanonicalForm& R, ///< [in,out] remainder, degree (R, 1) <
- ///< degree (G, 1)
- const CFList& MOD ///< [in] only contains powers of
- ///< Variables of level higher than 1
- );
-
-
-/// division with remainder of @a F by
-/// @a G wrt Variable (1) modulo @a M using Newton inversion
-///
-/// @return @a Q returns the dividend, @a R returns the remainder.
-/// @sa divrem2(), newtonDiv()
-void
-newtonDivrem (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
- const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
- ///< which is monic in Variable (1)
- CanonicalForm& Q, ///< [in,out] dividend
- CanonicalForm& R, ///< [in,out] remainder, degree (R, 1) <
- ///< degree (G, 1)
- const CanonicalForm& M ///< [in] power of Variable (2)
- );
-
-/// division of @a F by
-/// @a G wrt Variable (1) modulo @a M using Newton inversion
-///
-/// @return @a newtonDiv returns the dividend
-/// @sa divrem2(), newtonDivrem()
-CanonicalForm
-newtonDiv (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
- const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
- ///< which is monic in Variable (1)
- const CanonicalForm& M ///< [in] power of Variable (2)
- );
-
-/// reduce @a F modulo elements in @a M.
-///
-/// @return @a mod returns @a F modulo @a M
-CanonicalForm mod (const CanonicalForm& F, ///< [in] compressed polynomial
- const CFList& M ///< [in] list containing only
- ///< univariate polynomials
- );
-
-/// Karatsuba style modular multiplication for bivariate polynomials.
-///
-/// @return @a mulMod2 returns @a A * @a B mod @a M.
-CanonicalForm
-mulMod2 (const CanonicalForm& A, ///< [in] bivariate, compressed polynomial
- const CanonicalForm& B, ///< [in] bivariate, compressed polynomial
- const CanonicalForm& M ///< [in] power of Variable (2)
- );
-
-/// Karatsuba style modular multiplication for multivariate polynomials.
-///
-/// @return @a mulMod2 returns @a A * @a B mod @a MOD.
-CanonicalForm
-mulMod (const CanonicalForm& A, ///< [in] multivariate, compressed polynomial
- const CanonicalForm& B, ///< [in] multivariate, compressed polynomial
- const CFList& MOD ///< [in] only contains powers of
- ///< Variables of level higher than 1
- );
-
-/// product of all elements in @a L modulo @a M via divide-and-conquer.
-///
-/// @return @a prodMod returns product of all elements in @a L modulo @a M.
-CanonicalForm
-prodMod (const CFList& L, ///< [in] contains only bivariate, compressed
- ///< polynomials
- const CanonicalForm& M ///< [in] power of Variable (2)
- );
-
-/// product of all elements in @a L modulo @a M via divide-and-conquer.
-///
-/// @return @a prodMod returns product of all elements in @a L modulo @a M.
-CanonicalForm
-prodMod (const CFList& L, ///< [in] contains multivariate, compressed
- ///< polynomials
- const CFList& M ///< [in] contains only powers of Variables
- );
-
/// sort a list of polynomials by their degree in @a x.
///
void sortList (CFList& list, ///< [in, out] list of polys, sorted list
diff --git a/factory/facMul.cc b/factory/facMul.cc
new file mode 100644
index 0000000..f022ac5
--- /dev/null
+++ b/factory/facMul.cc
@@ -0,0 +1,2512 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR
+\*****************************************************************************/
+/** @file facHensel.cc
+ *
+ * This file implements functions for fast multiplication and division with
+ * remainder
+ *
+ * @author Martin Lee
+ *
+ **/
+/*****************************************************************************/
+
+#include "debug.h"
+
+#include "canonicalform.h"
+#include "facMul.h"
+#include "algext.h"
+#include "templates/ftmpl_functions.h"
+
+#ifdef HAVE_NTL
+#include <NTL/lzz_pEX.h>
+#include "NTLconvert.h"
+
+#ifdef HAVE_FLINT
+#include "FLINTconvert.h"
+#endif
+
+#ifdef HAVE_FLINT
+void kronSub (fmpz_poly_t result, const CanonicalForm& A, int d)
+{
+ int degAy= degree (A);
+ fmpz_poly_init2 (result, d*(degAy + 1));
+ _fmpz_poly_set_length (result, d*(degAy + 1));
+ CFIterator j;
+ for (CFIterator i= A; i.hasTerms(); i++)
+ {
+ if (i.coeff().inBaseDomain())
+ convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d), i.coeff());
+ else
+ for (j= i.coeff(); j.hasTerms(); j++)
+ convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d+j.exp()),
+ j.coeff());
+ }
+ _fmpz_poly_normalise(result);
+}
+
+
+CanonicalForm
+reverseSubst (const fmpz_poly_t F, int d, const Variable& alpha,
+ const CanonicalForm& den)
+{
+ Variable x= Variable (1);
+
+ CanonicalForm result= 0;
+ int i= 0;
+ int degf= fmpz_poly_degree (F);
+ int k= 0;
+ int degfSubK;
+ int repLength, j;
+ CanonicalForm coeff;
+ fmpz* tmp;
+ while (degf >= k)
+ {
+ coeff= 0;
+ degfSubK= degf - k;
+ if (degfSubK >= d)
+ repLength= d;
+ else
+ repLength= degfSubK + 1;
+
+ for (j= 0; j < repLength; j++)
+ {
+ tmp= fmpz_poly_get_coeff_ptr (F, j+k);
+ if (!fmpz_is_zero (tmp))
+ {
+ CanonicalForm ff= convertFmpz2CF (tmp)/den;
+ coeff += ff*power (alpha, j);
+ }
+ }
+ result += coeff*power (x, i);
+ i++;
+ k= d*i;
+ }
+ return result;
+}
+
+CanonicalForm
+mulFLINTQa (const CanonicalForm& F, const CanonicalForm& G,
+ const Variable& alpha)
+{
+ CanonicalForm A= F;
+ CanonicalForm B= G;
+
+ CanonicalForm denA= bCommonDen (A);
+ CanonicalForm denB= bCommonDen (B);
+
+ A *= denA;
+ B *= denB;
+ int degAa= degree (A, alpha);
+ int degBa= degree (B, alpha);
+ int d= degAa + 1 + degBa;
+
+ fmpz_poly_t FLINTA,FLINTB;
+ fmpz_poly_init (FLINTA);
+ fmpz_poly_init (FLINTB);
+ kronSub (FLINTA, A, d);
+ kronSub (FLINTB, B, d);
+
+ fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
+
+ denA *= denB;
+ A= reverseSubst (FLINTA, d, alpha, denA);
+
+ fmpz_poly_clear (FLINTA);
+ fmpz_poly_clear (FLINTB);
+ return A;
+}
+
+CanonicalForm
+mulFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
+{
+ CanonicalForm A= F;
+ CanonicalForm B= G;
+
+ CanonicalForm denA= bCommonDen (A);
+ CanonicalForm denB= bCommonDen (B);
+
+ A *= denA;
+ B *= denB;
+ fmpz_poly_t FLINTA,FLINTB;
+ convertFacCF2Fmpz_poly_t (FLINTA, A);
+ convertFacCF2Fmpz_poly_t (FLINTB, B);
+ fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
+ denA *= denB;
+ A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
+ A /= denA;
+ fmpz_poly_clear (FLINTA);
+ fmpz_poly_clear (FLINTB);
+
+ return A;
+}
+
+/*CanonicalForm
+mulFLINTQ2 (const CanonicalForm& F, const CanonicalForm& G)
+{
+ CanonicalForm A= F;
+ CanonicalForm B= G;
+
+ fmpq_poly_t FLINTA,FLINTB;
+ convertFacCF2Fmpq_poly_t (FLINTA, A);
+ convertFacCF2Fmpq_poly_t (FLINTB, B);
+
+ fmpq_poly_mul (FLINTA, FLINTA, FLINTB);
+ A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
+ fmpq_poly_clear (FLINTA);
+ fmpq_poly_clear (FLINTB);
+ return A;
+}*/
+
+CanonicalForm
+divFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
+{
+ CanonicalForm A= F;
+ CanonicalForm B= G;
+
+ fmpq_poly_t FLINTA,FLINTB;
+ convertFacCF2Fmpq_poly_t (FLINTA, A);
+ convertFacCF2Fmpq_poly_t (FLINTB, B);
+
+ fmpq_poly_div (FLINTA, FLINTA, FLINTB);
+ A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
+
+ fmpq_poly_clear (FLINTA);
+ fmpq_poly_clear (FLINTB);
+ return A;
+}
+
+CanonicalForm
+modFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
+{
+ CanonicalForm A= F;
+ CanonicalForm B= G;
+
+ fmpq_poly_t FLINTA,FLINTB;
+ convertFacCF2Fmpq_poly_t (FLINTA, A);
+ convertFacCF2Fmpq_poly_t (FLINTB, B);
+
+ fmpq_poly_rem (FLINTA, FLINTA, FLINTB);
+ A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
+
+ fmpq_poly_clear (FLINTA);
+ fmpq_poly_clear (FLINTB);
+ return A;
+}
+
+CanonicalForm
+mulFLINTQaTrunc (const CanonicalForm& F, const CanonicalForm& G,
+ const Variable& alpha, int m)
+{
+ CanonicalForm A= F;
+ CanonicalForm B= G;
+
+ CanonicalForm denA= bCommonDen (A);
+ CanonicalForm denB= bCommonDen (B);
+
+ A *= denA;
+ B *= denB;
+
+ int degAa= degree (A, alpha);
+ int degBa= degree (B, alpha);
+ int d= degAa + 1 + degBa;
+
+ fmpz_poly_t FLINTA,FLINTB;
+ fmpz_poly_init (FLINTA);
+ fmpz_poly_init (FLINTB);
+ kronSub (FLINTA, A, d);
+ kronSub (FLINTB, B, d);
+
+ int k= d*m;
+ fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, k);
+
+ denA *= denB;
+ A= reverseSubst (FLINTA, d, alpha, denA);
+ fmpz_poly_clear (FLINTA);
+ fmpz_poly_clear (FLINTB);
+ return A;
+}
+
+CanonicalForm
+mulFLINTQTrunc (const CanonicalForm& F, const CanonicalForm& G, int m)
+{
+ if (F.inCoeffDomain() || G.inCoeffDomain())
+ return mod (F*G, power (Variable (1), m));
+ Variable alpha;
+ if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
+ return mulFLINTQaTrunc (F, G, alpha, m);
+
+ CanonicalForm A= F;
+ CanonicalForm B= G;
+
+ CanonicalForm denA= bCommonDen (A);
+ CanonicalForm denB= bCommonDen (B);
+
+ A *= denA;
+ B *= denB;
+ fmpz_poly_t FLINTA,FLINTB;
+ convertFacCF2Fmpz_poly_t (FLINTA, A);
+ convertFacCF2Fmpz_poly_t (FLINTB, B);
+ fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, m);
+ denA *= denB;
+ A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
+ A /= denA;
+ fmpz_poly_clear (FLINTA);
+ fmpz_poly_clear (FLINTB);
+
+ return A;
+}
+
+CanonicalForm uniReverse (const CanonicalForm& F, int d)
+{
+ if (d == 0)
+ return F;
+ if (F.inCoeffDomain())
+ return F*power (Variable (1),d);
+ Variable x= Variable (1);
+ CanonicalForm result= 0;
+ CFIterator i= F;
+ while (d - i.exp() < 0)
+ i++;
+
+ for (; i.hasTerms() && (d - i.exp() >= 0); i++)
+ result += i.coeff()*power (x, d - i.exp());
+ return result;
+}
+
+CanonicalForm
+newtonInverse (const CanonicalForm& F, const int n)
+{
+ int l= ilog2(n);
+
+ CanonicalForm g= F [0];
+
+ ASSERT (!g.isZero(), "expected a unit");
+
+ if (!g.isOne())
+ g = 1/g;
+ Variable x= Variable (1);
+ CanonicalForm result;
+ int exp= 0;
+ if (n & 1)
+ {
+ result= g;
+ exp= 1;
+ }
+ CanonicalForm h;
+
+ for (int i= 1; i <= l; i++)
+ {
+ h= mulNTL (g, mod (F, power (x, (1 << i))));
+ h= mod (h, power (x, (1 << i)) - 1);
+ h= div (h, power (x, (1 << (i - 1))));
+ g -= power (x, (1 << (i - 1)))*
+ mulFLINTQTrunc (g, h, 1 << (i-1));
+
+ if (n & (1 << i))
+ {
+ if (exp)
+ {
+ h= mulNTL (result, mod (F, power (x, exp + (1 << i))));
+ h= mod (h, power (x, exp + (1 << i)) - 1);
+ h= div (h, power (x, exp));
+ result -= power(x, exp)*mulFLINTQTrunc (g, h, 1 << i);
+ exp += (1 << i);
+ }
+ else
+ {
+ exp= (1 << i);
+ result= g;
+ }
+ }
+ }
+
+ return result;
+}
+
+void
+newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+ CanonicalForm& R)
+{
+ CanonicalForm A= F;
+ CanonicalForm B= G;
+ Variable x= Variable (1);
+ int degA= degree (A, x);
+ int degB= degree (B, x);
+ int m= degA - degB;
+
+ if (m < 0)
+ {
+ R= A;
+ Q= 0;
+ return;
+ }
+
+ if (degB <= 1)
+ divrem (A, B, Q, R);
+ else
+ {
+ R= uniReverse (A, degA);
+
+ CanonicalForm revB= uniReverse (B, degB);
+ CanonicalForm buf= revB;
+ revB= newtonInverse (revB, m + 1);
+ Q= mulFLINTQTrunc (R, revB, m + 1);
+ Q= uniReverse (Q, m);
+
+ R= A - mulNTL (Q, B);
+ }
+}
+
+void
+newtonDiv (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q)
+{
+ CanonicalForm A= F;
+ CanonicalForm B= G;
+ Variable x= Variable (1);
+ int degA= degree (A, x);
+ int degB= degree (B, x);
+ int m= degA - degB;
+
+ if (m < 0)
+ {
+ Q= 0;
+ return;
+ }
+
+ if (degB <= 1)
+ Q= div (A, B);
+ else
+ {
+ CanonicalForm R= uniReverse (A, degA);
+
+ CanonicalForm revB= uniReverse (B, degB);
+ revB= newtonInverse (revB, m + 1);
+ Q= mulFLINTQTrunc (R, revB, m + 1);
+ Q= uniReverse (Q, m);
+ }
+}
+
+#endif
+
+CanonicalForm
+mulNTL (const CanonicalForm& F, const CanonicalForm& G)
+{
+ if (F.inCoeffDomain() || G.inCoeffDomain() || getCharacteristic() == 0)
+ {
+ Variable alpha;
+#ifdef HAVE_FLINT
+ if ((!F.inCoeffDomain() && !G.inCoeffDomain()) &&
+ (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha)))
+ {
+ CanonicalForm result= mulFLINTQa (F, G, alpha);
+ return result;
+ }
+ else if (!F.inCoeffDomain() && !G.inCoeffDomain())
+ return mulFLINTQ (F, G);
+#endif
+ return F*G;
+ }
+ ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
+ ASSERT (F.level() == G.level(), "expected polys of same level");
+ if (CFFactory::gettype() == GaloisFieldDomain)
+ return F*G;
+ zz_p::init (getCharacteristic());
+ Variable alpha;
+ CanonicalForm result;
+ if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
+ {
+ zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
+ zz_pE::init (NTLMipo);
+ zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
+ zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
+ mul (NTLF, NTLF, NTLG);
+ result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
+ }
+ else
+ {
+#ifdef HAVE_FLINT
+ nmod_poly_t FLINTF, FLINTG;
+ convertFacCF2nmod_poly_t (FLINTF, F);
+ convertFacCF2nmod_poly_t (FLINTG, G);
+ nmod_poly_mul (FLINTF, FLINTF, FLINTG);
+ result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
+ nmod_poly_clear (FLINTF);
+ nmod_poly_clear (FLINTG);
+#else
+ zz_pX NTLF= convertFacCF2NTLzzpX (F);
+ zz_pX NTLG= convertFacCF2NTLzzpX (G);
+ mul (NTLF, NTLF, NTLG);
+ result= convertNTLzzpX2CF(NTLF, F.mvar());
+#endif
+ }
+ return result;
+}
+
+CanonicalForm
+modNTL (const CanonicalForm& F, const CanonicalForm& G)
+{
+ if (F.inCoeffDomain() && G.isUnivariate())
+ return F;
+ else if (F.inCoeffDomain() && G.inCoeffDomain())
+ return mod (F, G);
+ else if (F.isUnivariate() && G.inCoeffDomain())
+ return mod (F,G);
+
+ if (getCharacteristic() == 0)
+ {
+#ifdef HAVE_FLINT
+ Variable alpha;
+ if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
+ return modFLINTQ (F, G);
+ else
+ {
+ CanonicalForm Q, R;
+ newtonDivrem (F, G, Q, R);
+ return R;
+ }
+#else
+ return mod (F, G);
+#endif
+ }
+
+ ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
+ ASSERT (F.level() == G.level(), "expected polys of same level");
+ if (CFFactory::gettype() == GaloisFieldDomain)
+ return mod (F, G);
+ zz_p::init (getCharacteristic());
+ Variable alpha;
+ CanonicalForm result;
+ if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
+ {
+ zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
+ zz_pE::init (NTLMipo);
+ zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
+ zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
+ rem (NTLF, NTLF, NTLG);
+ result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
+ }
+ else
+ {
+#ifdef HAVE_FLINT
+ nmod_poly_t FLINTF, FLINTG;
+ convertFacCF2nmod_poly_t (FLINTF, F);
+ convertFacCF2nmod_poly_t (FLINTG, G);
+ nmod_poly_divrem (FLINTG, FLINTF, FLINTF, FLINTG);
+ result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
+ nmod_poly_clear (FLINTF);
+ nmod_poly_clear (FLINTG);
+#else
+ zz_pX NTLF= convertFacCF2NTLzzpX (F);
+ zz_pX NTLG= convertFacCF2NTLzzpX (G);
+ rem (NTLF, NTLF, NTLG);
+ result= convertNTLzzpX2CF(NTLF, F.mvar());
+#endif
+ }
+ return result;
+}
+
+CanonicalForm
+divNTL (const CanonicalForm& F, const CanonicalForm& G)
+{
+ if (F.inCoeffDomain() && G.isUnivariate())
+ return F;
+ else if (F.inCoeffDomain() && G.inCoeffDomain())
+ return div (F, G);
+ else if (F.isUnivariate() && G.inCoeffDomain())
+ return div (F,G);
+
+ if (getCharacteristic() == 0)
+ {
+#ifdef HAVE_FLINT
+ Variable alpha;
+ if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
+ return divFLINTQ (F,G);
+ else
+ {
+ CanonicalForm Q;
+ newtonDiv (F, G, Q);
+ return Q;
+ }
+#else
+ return div (F, G);
+#endif
+ }
+
+ ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
+ ASSERT (F.level() == G.level(), "expected polys of same level");
+ if (CFFactory::gettype() == GaloisFieldDomain)
+ return div (F, G);
+ zz_p::init (getCharacteristic());
+ Variable alpha;
+ CanonicalForm result;
+ if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
+ {
+ zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
+ zz_pE::init (NTLMipo);
+ zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
+ zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
+ div (NTLF, NTLF, NTLG);
+ result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
+ }
+ else
+ {
+#ifdef HAVE_FLINT
+ nmod_poly_t FLINTF, FLINTG;
+ convertFacCF2nmod_poly_t (FLINTF, F);
+ convertFacCF2nmod_poly_t (FLINTG, G);
+ nmod_poly_div (FLINTF, FLINTF, FLINTG);
+ result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
+ nmod_poly_clear (FLINTF);
+ nmod_poly_clear (FLINTG);
+#else
+ zz_pX NTLF= convertFacCF2NTLzzpX (F);
+ zz_pX NTLG= convertFacCF2NTLzzpX (G);
+ div (NTLF, NTLF, NTLG);
+ result= convertNTLzzpX2CF(NTLF, F.mvar());
+#endif
+ }
+ return result;
+}
+
+#ifdef HAVE_FLINT
+void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
+{
+ int degAy= degree (A);
+ nmod_poly_init2 (result, getCharacteristic(), d*(degAy + 1));
+
+ nmod_poly_t buf;
+
+ int j, k, bufRepLength;
+ for (CFIterator i= A; i.hasTerms(); i++)
+ {
+ convertFacCF2nmod_poly_t (buf, i.coeff());
+
+ k= i.exp()*d;
+ bufRepLength= (int) nmod_poly_length (buf);
+ for (j= 0; j < bufRepLength; j++)
+ nmod_poly_set_coeff_ui (result, j + k, nmod_poly_get_coeff_ui (buf, j));
+ nmod_poly_clear (buf);
+ }
+ _nmod_poly_normalise (result);
+}
+
+/*void kronSubQ (fmpz_poly_t result, const CanonicalForm& A, int d)
+{
+ int degAy= degree (A);
+ fmpz_poly_init2 (result, d*(degAy + 1));
+ _fmpz_poly_set_length (result, d*(degAy+1));
+
+ CFIterator j;
+ for (CFIterator i= A; i.hasTerms(); i++)
+ {
+ if (i.coeff().inBas
+ convertFacCF2Fmpz_poly_t (buf, i.coeff());
+
+ int k= i.exp()*d;
+ int bufRepLength= (int) fmpz_poly_length (buf);
+ for (int j= 0; j < bufRepLength; j++)
+ {
+ fmpz_poly_get_coeff_fmpz (coeff, buf, j);
+ fmpz_poly_set_coeff_fmpz (result, j + k, coeff);
+ }
+ fmpz_poly_clear (buf);
+ }
+ fmpz_clear (coeff);
+ _fmpz_poly_normalise (result);
+}*/
+
+// A is a bivariate poly over Qa!!!!
+// d2= 2*deg_alpha + 1
+// d1= 2*deg_x*d2+1
+void kronSubQa (fmpq_poly_t result, const CanonicalForm& A, int d1, int d2)
+{
+ int degAy= degree (A);
+ fmpq_poly_init2 (result, d1*(degAy + 1));
+
+ fmpq_poly_t buf;
+ fmpq_t coeff;
+
+ int k, l, bufRepLength;
+ CFIterator j;
+ for (CFIterator i= A; i.hasTerms(); i++)
+ {
+ if (i.coeff().inCoeffDomain())
+ {
+ k= d1*i.exp();
+ convertFacCF2Fmpq_poly_t (buf, i.coeff());
+ bufRepLength= (int) fmpq_poly_length(buf);
+ for (l= 0; l < bufRepLength; l++)
+ {
+ fmpq_poly_get_coeff_fmpq (coeff, buf, l);
+ fmpq_poly_set_coeff_fmpq (result, l + k, coeff);
+ }
+ fmpq_poly_clear (buf);
+ }
+ else
+ {
+ for (j= i.coeff(); j.hasTerms(); j++)
+ {
+ k= d1*i.exp();
+ k += d2*j.exp();
+ convertFacCF2Fmpq_poly_t (buf, j.coeff());
+ bufRepLength= (int) fmpq_poly_length(buf);
+ for (l= 0; l < bufRepLength; l++)
+ {
+ fmpq_poly_get_coeff_fmpq (coeff, buf, l);
+ fmpq_poly_set_coeff_fmpq (result, k + l, coeff);
+ }
+ fmpq_poly_clear (buf);
+ }
+ }
+ }
+ fmpq_clear (coeff);
+ _fmpq_poly_normalise (result);
+}
+#endif
+
+zz_pX kronSubFp (const CanonicalForm& A, int d)
+{
+ int degAy= degree (A);
+ zz_pX result;
+ result.rep.SetLength (d*(degAy + 1));
+
+ zz_p *resultp;
+ resultp= result.rep.elts();
+ zz_pX buf;
+ zz_p *bufp;
+ int j, k, bufRepLength;
+
+ for (CFIterator i= A; i.hasTerms(); i++)
+ {
+ if (i.coeff().inCoeffDomain())
+ buf= convertFacCF2NTLzzpX (i.coeff());
+ else
+ buf= convertFacCF2NTLzzpX (i.coeff());
+
+ k= i.exp()*d;
+ bufp= buf.rep.elts();
+ bufRepLength= (int) buf.rep.length();
+ for (j= 0; j < bufRepLength; j++)
+ resultp [j + k]= bufp [j];
+ }
+ result.normalize();
+
+ return result;
+}
+
+zz_pEX kronSub (const CanonicalForm& A, int d, const Variable& alpha)
+{
+ int degAy= degree (A);
+ zz_pEX result;
+ result.rep.SetLength (d*(degAy + 1));
+
+ Variable v;
+ zz_pE *resultp;
+ resultp= result.rep.elts();
+ zz_pEX buf1;
+ zz_pE *buf1p;
+ zz_pX buf2;
+ zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
+ int j, k, buf1RepLength;
+
+ for (CFIterator i= A; i.hasTerms(); i++)
+ {
+ if (i.coeff().inCoeffDomain())
+ {
+ buf2= convertFacCF2NTLzzpX (i.coeff());
+ buf1= to_zz_pEX (to_zz_pE (buf2));
+ }
+ else
+ buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
+
+ k= i.exp()*d;
+ buf1p= buf1.rep.elts();
+ buf1RepLength= (int) buf1.rep.length();
+ for (j= 0; j < buf1RepLength; j++)
+ resultp [j + k]= buf1p [j];
+ }
+ result.normalize();
+
+ return result;
+}
+
+void
+kronSubRecipro (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
+ const Variable& alpha)
+{
+ int degAy= degree (A);
+ subA1.rep.SetLength ((long) d*(degAy + 2));
+ subA2.rep.SetLength ((long) d*(degAy + 2));
+
+ Variable v;
+ zz_pE *subA1p;
+ zz_pE *subA2p;
+ subA1p= subA1.rep.elts();
+ subA2p= subA2.rep.elts();
+ zz_pEX buf;
+ zz_pE *bufp;
+ zz_pX buf2;
+ zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
+ int j, k, kk, bufRepLength;
+
+ for (CFIterator i= A; i.hasTerms(); i++)
+ {
+ if (i.coeff().inCoeffDomain())
+ {
+ buf2= convertFacCF2NTLzzpX (i.coeff());
+ buf= to_zz_pEX (to_zz_pE (buf2));
+ }
+ else
+ buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
+
+ k= i.exp()*d;
+ kk= (degAy - i.exp())*d;
+ bufp= buf.rep.elts();
+ bufRepLength= (int) buf.rep.length();
+ for (j= 0; j < bufRepLength; j++)
+ {
+ subA1p [j + k] += bufp [j];
+ subA2p [j + kk] += bufp [j];
+ }
+ }
+ subA1.normalize();
+ subA2.normalize();
+}
+
+void
+kronSubRecipro (zz_pX& subA1, zz_pX& subA2,const CanonicalForm& A, int d)
+{
+ int degAy= degree (A);
+ subA1.rep.SetLength ((long) d*(degAy + 2));
+ subA2.rep.SetLength ((long) d*(degAy + 2));
+
+ zz_p *subA1p;
+ zz_p *subA2p;
+ subA1p= subA1.rep.elts();
+ subA2p= subA2.rep.elts();
+ zz_pX buf;
+ zz_p *bufp;
+ int j, k, kk, bufRepLength;
+
+ for (CFIterator i= A; i.hasTerms(); i++)
+ {
+ buf= convertFacCF2NTLzzpX (i.coeff());
+
+ k= i.exp()*d;
+ kk= (degAy - i.exp())*d;
+ bufp= buf.rep.elts();
+ bufRepLength= (int) buf.rep.length();
+ for (j= 0; j < bufRepLength; j++)
+ {
+ subA1p [j + k] += bufp [j];
+ subA2p [j + kk] += bufp [j];
+ }
+ }
+ subA1.normalize();
+ subA2.normalize();
+}
+
+CanonicalForm
+reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
+ const Variable& alpha)
+{
+ Variable y= Variable (2);
+ Variable x= Variable (1);
+
+ zz_pEX f= F;
+ zz_pEX g= G;
+ int degf= deg(f);
+ int degg= deg(g);
+
+ zz_pEX buf1;
+ zz_pEX buf2;
+ zz_pEX buf3;
+
+ zz_pE *buf1p;
+ zz_pE *buf2p;
+ zz_pE *buf3p;
+ if (f.rep.length() < (long) d*(k+1)) //zero padding
+ f.rep.SetLength ((long)d*(k+1));
+
+ zz_pE *gp= g.rep.elts();
+ zz_pE *fp= f.rep.elts();
+ CanonicalForm result= 0;
+ int i= 0;
+ int lf= 0;
+ int lg= d*k;
+ int degfSubLf= degf;
+ int deggSubLg= degg-lg;
+ int repLengthBuf2, repLengthBuf1, ind, tmp;
+ zz_pE zzpEZero= zz_pE();
+
+ while (degf >= lf || lg >= 0)
+ {
+ if (degfSubLf >= d)
+ repLengthBuf1= d;
+ else if (degfSubLf < 0)
+ repLengthBuf1= 0;
+ else
+ repLengthBuf1= degfSubLf + 1;
+ buf1.rep.SetLength((long) repLengthBuf1);
+
+ buf1p= buf1.rep.elts();
+ for (ind= 0; ind < repLengthBuf1; ind++)
+ buf1p [ind]= fp [ind + lf];
+ buf1.normalize();
+
+ repLengthBuf1= buf1.rep.length();
+
+ if (deggSubLg >= d - 1)
+ repLengthBuf2= d - 1;
+ else if (deggSubLg < 0)
+ repLengthBuf2= 0;
+ else
+ repLengthBuf2= deggSubLg + 1;
+
+ buf2.rep.SetLength ((long) repLengthBuf2);
+ buf2p= buf2.rep.elts();
+ for (ind= 0; ind < repLengthBuf2; ind++)
+ buf2p [ind]= gp [ind + lg];
+ buf2.normalize();
+
+ repLengthBuf2= buf2.rep.length();
+
+ buf3.rep.SetLength((long) repLengthBuf2 + d);
+ buf3p= buf3.rep.elts();
+ buf2p= buf2.rep.elts();
+ buf1p= buf1.rep.elts();
+ for (ind= 0; ind < repLengthBuf1; ind++)
+ buf3p [ind]= buf1p [ind];
+ for (ind= repLengthBuf1; ind < d; ind++)
+ buf3p [ind]= zzpEZero;
+ for (ind= 0; ind < repLengthBuf2; ind++)
+ buf3p [ind + d]= buf2p [ind];
+ buf3.normalize();
+
+ result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
+ i++;
+
+
+ lf= i*d;
+ degfSubLf= degf - lf;
+
+ lg= d*(k-i);
+ deggSubLg= degg - lg;
+
+ buf1p= buf1.rep.elts();
+
+ if (lg >= 0 && deggSubLg > 0)
+ {
+ if (repLengthBuf2 > degfSubLf + 1)
+ degfSubLf= repLengthBuf2 - 1;
+ tmp= tmin (repLengthBuf1, deggSubLg + 1);
+ for (ind= 0; ind < tmp; ind++)
+ gp [ind + lg] -= buf1p [ind];
+ }
+
+ if (lg < 0)
+ break;
+
+ buf2p= buf2.rep.elts();
+ if (degfSubLf >= 0)
+ {
+ for (ind= 0; ind < repLengthBuf2; ind++)
+ fp [ind + lf] -= buf2p [ind];
+ }
+ }
+
+ return result;
+}
+
+CanonicalForm
+reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
+{
+ Variable y= Variable (2);
+ Variable x= Variable (1);
+
+ zz_pX f= F;
+ zz_pX g= G;
+ int degf= deg(f);
+ int degg= deg(g);
+
+ zz_pX buf1;
+ zz_pX buf2;
+ zz_pX buf3;
+
+ zz_p *buf1p;
+ zz_p *buf2p;
+ zz_p *buf3p;
+
+ if (f.rep.length() < (long) d*(k+1)) //zero padding
+ f.rep.SetLength ((long)d*(k+1));
+
+ zz_p *gp= g.rep.elts();
+ zz_p *fp= f.rep.elts();
+ CanonicalForm result= 0;
+ int i= 0;
+ int lf= 0;
+ int lg= d*k;
+ int degfSubLf= degf;
+ int deggSubLg= degg-lg;
+ int repLengthBuf2, repLengthBuf1, ind, tmp;
+ zz_p zzpZero= zz_p();
+ while (degf >= lf || lg >= 0)
+ {
+ if (degfSubLf >= d)
+ repLengthBuf1= d;
+ else if (degfSubLf < 0)
+ repLengthBuf1= 0;
+ else
+ repLengthBuf1= degfSubLf + 1;
+ buf1.rep.SetLength((long) repLengthBuf1);
+
+ buf1p= buf1.rep.elts();
+ for (ind= 0; ind < repLengthBuf1; ind++)
+ buf1p [ind]= fp [ind + lf];
+ buf1.normalize();
+
+ repLengthBuf1= buf1.rep.length();
+
+ if (deggSubLg >= d - 1)
+ repLengthBuf2= d - 1;
+ else if (deggSubLg < 0)
+ repLengthBuf2= 0;
+ else
+ repLengthBuf2= deggSubLg + 1;
+
+ buf2.rep.SetLength ((long) repLengthBuf2);
+ buf2p= buf2.rep.elts();
+ for (ind= 0; ind < repLengthBuf2; ind++)
+ buf2p [ind]= gp [ind + lg];
+
+ buf2.normalize();
+
+ repLengthBuf2= buf2.rep.length();
+
+
+ buf3.rep.SetLength((long) repLengthBuf2 + d);
+ buf3p= buf3.rep.elts();
+ buf2p= buf2.rep.elts();
+ buf1p= buf1.rep.elts();
+ for (ind= 0; ind < repLengthBuf1; ind++)
+ buf3p [ind]= buf1p [ind];
+ for (ind= repLengthBuf1; ind < d; ind++)
+ buf3p [ind]= zzpZero;
+ for (ind= 0; ind < repLengthBuf2; ind++)
+ buf3p [ind + d]= buf2p [ind];
+ buf3.normalize();
+
+ result += convertNTLzzpX2CF (buf3, x)*power (y, i);
+ i++;
+
+
+ lf= i*d;
+ degfSubLf= degf - lf;
+
+ lg= d*(k-i);
+ deggSubLg= degg - lg;
+
+ buf1p= buf1.rep.elts();
+
+ if (lg >= 0 && deggSubLg > 0)
+ {
+ if (repLengthBuf2 > degfSubLf + 1)
+ degfSubLf= repLengthBuf2 - 1;
+ tmp= tmin (repLengthBuf1, deggSubLg + 1);
+ for (ind= 0; ind < tmp; ind++)
+ gp [ind + lg] -= buf1p [ind];
+ }
+ if (lg < 0)
+ break;
+
+ buf2p= buf2.rep.elts();
+ if (degfSubLf >= 0)
+ {
+ for (ind= 0; ind < repLengthBuf2; ind++)
+ fp [ind + lf] -= buf2p [ind];
+ }
+ }
+
+ return result;
+}
+
+CanonicalForm reverseSubst (const zz_pEX& F, int d, const Variable& alpha)
+{
+ Variable y= Variable (2);
+ Variable x= Variable (1);
+
+ zz_pEX f= F;
+ zz_pE *fp= f.rep.elts();
+
+ zz_pEX buf;
+ zz_pE *bufp;
+ CanonicalForm result= 0;
+ int i= 0;
+ int degf= deg(f);
+ int k= 0;
+ int degfSubK, repLength, j;
+ while (degf >= k)
+ {
+ degfSubK= degf - k;
+ if (degfSubK >= d)
+ repLength= d;
+ else
+ repLength= degfSubK + 1;
+
+ buf.rep.SetLength ((long) repLength);
+ bufp= buf.rep.elts();
+ for (j= 0; j < repLength; j++)
+ bufp [j]= fp [j + k];
+ buf.normalize();
+
+ result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
+ i++;
+ k= d*i;
+ }
+
+ return result;
+}
+
+CanonicalForm reverseSubstFp (const zz_pX& F, int d)
+{
+ Variable y= Variable (2);
+ Variable x= Variable (1);
+
+ zz_pX f= F;
+ zz_p *fp= f.rep.elts();
+
+ zz_pX buf;
+ zz_p *bufp;
+ CanonicalForm result= 0;
+ int i= 0;
+ int degf= deg(f);
+ int k= 0;
+ int degfSubK, repLength, j;
+ while (degf >= k)
+ {
+ degfSubK= degf - k;
+ if (degfSubK >= d)
+ repLength= d;
+ else
+ repLength= degfSubK + 1;
+
+ buf.rep.SetLength ((long) repLength);
+ bufp= buf.rep.elts();
+ for (j= 0; j < repLength; j++)
+ bufp [j]= fp [j + k];
+ buf.normalize();
+
+ result += convertNTLzzpX2CF (buf, x)*power (y, i);
+ i++;
+ k= d*i;
+ }
+
+ return result;
+}
+
+// assumes input to be reduced mod M and to be an element of Fq not Fp
+CanonicalForm
+mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
+ CanonicalForm& M)
+{
+ int d1= degree (F, 1) + degree (G, 1) + 1;
+ d1 /= 2;
+ d1 += 1;
+
+ zz_pX F1, F2;
+ kronSubRecipro (F1, F2, F, d1);
+ zz_pX G1, G2;
+ kronSubRecipro (G1, G2, G, d1);
+
+ int k= d1*degree (M);
+ MulTrunc (F1, F1, G1, (long) k);
+
+ int degtailF= degree (tailcoeff (F), 1);
+ int degtailG= degree (tailcoeff (G), 1);
+ int taildegF= taildegree (F);
+ int taildegG= taildegree (G);
+ int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
+
+ reverse (F2, F2);
+ reverse (G2, G2);
+ MulTrunc (F2, F2, G2, b + 1);
+ reverse (F2, F2, b);
+
+ int d2= tmax (deg (F2)/d1, deg (F1)/d1);
+ return reverseSubst (F1, F2, d1, d2);
+}
+
+//Kronecker substitution
+CanonicalForm
+mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
+ CanonicalForm& M)
+{
+ CanonicalForm A= F;
+ CanonicalForm B= G;
+
+ int degAx= degree (A, 1);
+ int degAy= degree (A, 2);
+ int degBx= degree (B, 1);
+ int degBy= degree (B, 2);
+ int d1= degAx + 1 + degBx;
+ int d2= tmax (degAy, degBy);
+
+ if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
+ return mulMod2NTLFpReci (A, B, M);
+
+ zz_pX NTLA= kronSubFp (A, d1);
+ zz_pX NTLB= kronSubFp (B, d1);
+
+ int k= d1*degree (M);
+ MulTrunc (NTLA, NTLA, NTLB, (long) k);
+
+ A= reverseSubstFp (NTLA, d1);
+
+ return A;
+}
+
+// assumes input to be reduced mod M and to be an element of Fq not Fp
+CanonicalForm
+mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
+ CanonicalForm& M, const Variable& alpha)
+{
+ int d1= degree (F, 1) + degree (G, 1) + 1;
+ d1 /= 2;
+ d1 += 1;
+
+ zz_pEX F1, F2;
+ kronSubRecipro (F1, F2, F, d1, alpha);
+ zz_pEX G1, G2;
+ kronSubRecipro (G1, G2, G, d1, alpha);
+
+ int k= d1*degree (M);
+ MulTrunc (F1, F1, G1, (long) k);
+
+ int degtailF= degree (tailcoeff (F), 1);
+ int degtailG= degree (tailcoeff (G), 1);
+ int taildegF= taildegree (F);
+ int taildegG= taildegree (G);
+ int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
+
+ reverse (F2, F2);
+ reverse (G2, G2);
+ MulTrunc (F2, F2, G2, b + 1);
+ reverse (F2, F2, b);
+
+ int d2= tmax (deg (F2)/d1, deg (F1)/d1);
+ return reverseSubst (F1, F2, d1, d2, alpha);
+}
+
+#ifdef HAVE_FLINT
+CanonicalForm
+mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
+ CanonicalForm& M);
+#endif
+
+CanonicalForm
+mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
+ CanonicalForm& M)
+{
+ Variable alpha;
+ CanonicalForm A= F;
+ CanonicalForm B= G;
+
+ if (hasFirstAlgVar (A, alpha) || hasFirstAlgVar (B, alpha))
+ {
+ int degAx= degree (A, 1);
+ int degAy= degree (A, 2);
+ int degBx= degree (B, 1);
+ int degBy= degree (B, 2);
+ int d1= degAx + degBx + 1;
+ int d2= tmax (degAy, degBy);
+ zz_p::init (getCharacteristic());
+ zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
+ zz_pE::init (NTLMipo);
+
+ int degMipo= degree (getMipo (alpha));
+ if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy) &&
+ (2*degAy > degree (M)))
+ return mulMod2NTLFqReci (A, B, M, alpha);
+
+ zz_pEX NTLA= kronSub (A, d1, alpha);
+ zz_pEX NTLB= kronSub (B, d1, alpha);
+
+ int k= d1*degree (M);
+
+ MulTrunc (NTLA, NTLA, NTLB, (long) k);
+
+ A= reverseSubst (NTLA, d1, alpha);
+
+ return A;
+ }
+ else
+#ifdef HAVE_FLINT
+ return mulMod2FLINTFp (A, B, M);
+#else
+ return mulMod2NTLFp (A, B, M);
+#endif
+}
+
+#ifdef HAVE_FLINT
+void
+kronSubRecipro (nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm& A,
+ int d)
+{
+ int degAy= degree (A);
+ mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
+ nmod_poly_init2_preinv (subA1, getCharacteristic(), ninv, d*(degAy + 2));
+ nmod_poly_init2_preinv (subA2, getCharacteristic(), ninv, d*(degAy + 2));
+
+ nmod_poly_t buf;
+
+ int k, kk, j, bufRepLength;
+ for (CFIterator i= A; i.hasTerms(); i++)
+ {
+ convertFacCF2nmod_poly_t (buf, i.coeff());
+
+ k= i.exp()*d;
+ kk= (degAy - i.exp())*d;
+ bufRepLength= (int) nmod_poly_length (buf);
+ for (j= 0; j < bufRepLength; j++)
+ {
+ nmod_poly_set_coeff_ui (subA1, j + k,
+ n_addmod (nmod_poly_get_coeff_ui (subA1, j+k),
+ nmod_poly_get_coeff_ui (buf, j),
+ getCharacteristic()
+ )
+ );
+ nmod_poly_set_coeff_ui (subA2, j + kk,
+ n_addmod (nmod_poly_get_coeff_ui (subA2, j + kk),
+ nmod_poly_get_coeff_ui (buf, j),
+ getCharacteristic()
+ )
+ );
+ }
+ nmod_poly_clear (buf);
+ }
+ _nmod_poly_normalise (subA1);
+ _nmod_poly_normalise (subA2);
+}
+
+void
+kronSubRecipro (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
+ int d)
+{
+ int degAy= degree (A);
+ fmpz_poly_init2 (subA1, d*(degAy + 2));
+ fmpz_poly_init2 (subA2, d*(degAy + 2));
+
+ fmpz_poly_t buf;
+ fmpz_t coeff1, coeff2;
+
+ int k, kk, j, bufRepLength;
+ for (CFIterator i= A; i.hasTerms(); i++)
+ {
+ convertFacCF2Fmpz_poly_t (buf, i.coeff());
+
+ k= i.exp()*d;
+ kk= (degAy - i.exp())*d;
+ bufRepLength= (int) fmpz_poly_length (buf);
+ for (j= 0; j < bufRepLength; j++)
+ {
+ fmpz_poly_get_coeff_fmpz (coeff1, subA1, j+k);
+ fmpz_poly_get_coeff_fmpz (coeff2, buf, j);
+ fmpz_add (coeff1, coeff1, coeff2);
+ fmpz_poly_set_coeff_fmpz (subA1, j + k, coeff1);
+ fmpz_poly_get_coeff_fmpz (coeff1, subA2, j + kk);
+ fmpz_add (coeff1, coeff1, coeff2);
+ fmpz_poly_set_coeff_fmpz (subA2, j + kk, coeff1);
+ }
+ fmpz_poly_clear (buf);
+ }
+ fmpz_clear (coeff1);
+ fmpz_clear (coeff2);
+ _fmpz_poly_normalise (subA1);
+ _fmpz_poly_normalise (subA2);
+}
+
+CanonicalForm reverseSubstQ (const fmpz_poly_t F, int d)
+{
+ Variable y= Variable (2);
+ Variable x= Variable (1);
+
+ fmpz_poly_t f;
+ fmpz_poly_init (f);
+ fmpz_poly_set (f, F);
+
+ fmpz_poly_t buf;
+ CanonicalForm result= 0;
+ int i= 0;
+ int degf= fmpz_poly_degree(f);
+ int k= 0;
+ int degfSubK, repLength, j;
+ fmpz_t coeff;
+ while (degf >= k)
+ {
+ degfSubK= degf - k;
+ if (degfSubK >= d)
+ repLength= d;
+ else
+ repLength= degfSubK + 1;
+
+ fmpz_poly_init2 (buf, repLength);
+ fmpz_init (coeff);
+ for (j= 0; j < repLength; j++)
+ {
+ fmpz_poly_get_coeff_fmpz (coeff, f, j + k);
+ fmpz_poly_set_coeff_fmpz (buf, j, coeff);
+ }
+ _fmpz_poly_normalise (buf);
+
+ result += convertFmpz_poly_t2FacCF (buf, x)*power (y, i);
+ i++;
+ k= d*i;
+ fmpz_poly_clear (buf);
+ fmpz_clear (coeff);
+ }
+ fmpz_poly_clear (f);
+
+ return result;
+}
+
+CanonicalForm
+reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
+{
+ Variable y= Variable (2);
+ Variable x= Variable (1);
+
+ nmod_poly_t f, g;
+ mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
+ nmod_poly_init_preinv (f, getCharacteristic(), ninv);
+ nmod_poly_init_preinv (g, getCharacteristic(), ninv);
+ nmod_poly_set (f, F);
+ nmod_poly_set (g, G);
+ int degf= nmod_poly_degree(f);
+ int degg= nmod_poly_degree(g);
+
+
+ nmod_poly_t buf1,buf2, buf3;
+
+ if (nmod_poly_length (f) < (long) d*(k+1)) //zero padding
+ nmod_poly_fit_length (f,(long)d*(k+1));
+
+ CanonicalForm result= 0;
+ int i= 0;
+ int lf= 0;
+ int lg= d*k;
+ int degfSubLf= degf;
+ int deggSubLg= degg-lg;
+ int repLengthBuf2, repLengthBuf1, ind, tmp;
+ while (degf >= lf || lg >= 0)
+ {
+ if (degfSubLf >= d)
+ repLengthBuf1= d;
+ else if (degfSubLf < 0)
+ repLengthBuf1= 0;
+ else
+ repLengthBuf1= degfSubLf + 1;
+ nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
+
+ for (ind= 0; ind < repLengthBuf1; ind++)
+ nmod_poly_set_coeff_ui (buf1, ind, nmod_poly_get_coeff_ui (f, ind+lf));
+ _nmod_poly_normalise (buf1);
+
+ repLengthBuf1= nmod_poly_length (buf1);
+
+ if (deggSubLg >= d - 1)
+ repLengthBuf2= d - 1;
+ else if (deggSubLg < 0)
+ repLengthBuf2= 0;
+ else
+ repLengthBuf2= deggSubLg + 1;
+
+ nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
+ for (ind= 0; ind < repLengthBuf2; ind++)
+ nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
+
+ _nmod_poly_normalise (buf2);
+ repLengthBuf2= nmod_poly_length (buf2);
+
+ nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
+ for (ind= 0; ind < repLengthBuf1; ind++)
+ nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
+ for (ind= repLengthBuf1; ind < d; ind++)
+ nmod_poly_set_coeff_ui (buf3, ind, 0);
+ for (ind= 0; ind < repLengthBuf2; ind++)
+ nmod_poly_set_coeff_ui (buf3, ind+d, nmod_poly_get_coeff_ui (buf2, ind));
+ _nmod_poly_normalise (buf3);
+
+ result += convertnmod_poly_t2FacCF (buf3, x)*power (y, i);
+ i++;
+
+
+ lf= i*d;
+ degfSubLf= degf - lf;
+
+ lg= d*(k-i);
+ deggSubLg= degg - lg;
+
+ if (lg >= 0 && deggSubLg > 0)
+ {
+ if (repLengthBuf2 > degfSubLf + 1)
+ degfSubLf= repLengthBuf2 - 1;
+ tmp= tmin (repLengthBuf1, deggSubLg + 1);
+ for (ind= 0; ind < tmp; ind++)
+ nmod_poly_set_coeff_ui (g, ind + lg,
+ n_submod (nmod_poly_get_coeff_ui (g, ind + lg),
+ nmod_poly_get_coeff_ui (buf1, ind),
+ getCharacteristic()
+ )
+ );
+ }
+ if (lg < 0)
+ {
+ nmod_poly_clear (buf1);
+ nmod_poly_clear (buf2);
+ nmod_poly_clear (buf3);
+ break;
+ }
+ if (degfSubLf >= 0)
+ {
+ for (ind= 0; ind < repLengthBuf2; ind++)
+ nmod_poly_set_coeff_ui (f, ind + lf,
+ n_submod (nmod_poly_get_coeff_ui (f, ind + lf),
+ nmod_poly_get_coeff_ui (buf2, ind),
+ getCharacteristic()
+ )
+ );
+ }
+ nmod_poly_clear (buf1);
+ nmod_poly_clear (buf2);
+ nmod_poly_clear (buf3);
+ }
+
+ nmod_poly_clear (f);
+ nmod_poly_clear (g);
+
+ return result;
+}
+
+CanonicalForm
+reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
+{
+ Variable y= Variable (2);
+ Variable x= Variable (1);
+
+ fmpz_poly_t f, g;
+ fmpz_poly_init (f);
+ fmpz_poly_init (g);
+ fmpz_poly_set (f, F);
+ fmpz_poly_set (g, G);
+ int degf= fmpz_poly_degree(f);
+ int degg= fmpz_poly_degree(g);
+
+
+ fmpz_poly_t buf1,buf2, buf3;
+
+ if (fmpz_poly_length (f) < (long) d*(k+1)) //zero padding
+ fmpz_poly_fit_length (f,(long)d*(k+1));
+
+ CanonicalForm result= 0;
+ int i= 0;
+ int lf= 0;
+ int lg= d*k;
+ int degfSubLf= degf;
+ int deggSubLg= degg-lg;
+ int repLengthBuf2, repLengthBuf1, ind, tmp;
+ fmpz_t tmp1, tmp2;
+ while (degf >= lf || lg >= 0)
+ {
+ if (degfSubLf >= d)
+ repLengthBuf1= d;
+ else if (degfSubLf < 0)
+ repLengthBuf1= 0;
+ else
+ repLengthBuf1= degfSubLf + 1;
+
+ fmpz_poly_init2 (buf1, repLengthBuf1);
+
+ for (ind= 0; ind < repLengthBuf1; ind++)
+ {
+ fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
+ fmpz_poly_set_coeff_fmpz (buf1, ind, tmp1);
+ }
+ _fmpz_poly_normalise (buf1);
+
+ repLengthBuf1= fmpz_poly_length (buf1);
+
+ if (deggSubLg >= d - 1)
+ repLengthBuf2= d - 1;
+ else if (deggSubLg < 0)
+ repLengthBuf2= 0;
+ else
+ repLengthBuf2= deggSubLg + 1;
+
+ fmpz_poly_init2 (buf2, repLengthBuf2);
+
+ for (ind= 0; ind < repLengthBuf2; ind++)
+ {
+ fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
+ fmpz_poly_set_coeff_fmpz (buf2, ind, tmp1);
+ }
+
+ _fmpz_poly_normalise (buf2);
+ repLengthBuf2= fmpz_poly_length (buf2);
+
+ fmpz_poly_init2 (buf3, repLengthBuf2 + d);
+ for (ind= 0; ind < repLengthBuf1; ind++)
+ {
+ 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))
+ fmpz_poly_set_coeff_fmpz (buf3, ind, tmp1);
+ }
+ for (ind= repLengthBuf1; ind < d; ind++)
+ fmpz_poly_set_coeff_ui (buf3, ind, 0);
+ for (ind= 0; ind < repLengthBuf2; ind++)
+ {
+ fmpz_poly_get_coeff_fmpz (tmp1, buf2, ind);
+ fmpz_poly_set_coeff_fmpz (buf3, ind + d, tmp1);
+ }
+ _fmpz_poly_normalise (buf3);
+
+ result += convertFmpz_poly_t2FacCF (buf3, x)*power (y, i);
+ i++;
+
+
+ lf= i*d;
+ degfSubLf= degf - lf;
+
+ lg= d*(k-i);
+ deggSubLg= degg - lg;
+
+ if (lg >= 0 && deggSubLg > 0)
+ {
+ if (repLengthBuf2 > degfSubLf + 1)
+ degfSubLf= repLengthBuf2 - 1;
+ tmp= tmin (repLengthBuf1, deggSubLg + 1);
+ for (ind= 0; ind < tmp; ind++)
+ {
+ fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
+ fmpz_poly_get_coeff_fmpz (tmp2, buf1, ind);
+ fmpz_sub (tmp1, tmp1, tmp2);
+ fmpz_poly_set_coeff_fmpz (g, ind + lg, tmp1);
+ }
+ }
+ if (lg < 0)
+ {
+ fmpz_poly_clear (buf1);
+ fmpz_poly_clear (buf2);
+ fmpz_poly_clear (buf3);
+ break;
+ }
+ if (degfSubLf >= 0)
+ {
+ for (ind= 0; ind < repLengthBuf2; ind++)
+ {
+ fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
+ fmpz_poly_get_coeff_fmpz (tmp2, buf2, ind);
+ fmpz_sub (tmp1, tmp1, tmp2);
+ fmpz_poly_set_coeff_fmpz (f, ind + lf, tmp1);
+ }
+ }
+ fmpz_poly_clear (buf1);
+ fmpz_poly_clear (buf2);
+ fmpz_poly_clear (buf3);
+ }
+
+ fmpz_poly_clear (f);
+ fmpz_poly_clear (g);
+ fmpz_clear (tmp1);
+ fmpz_clear (tmp2);
+
+ return result;
+}
+
+CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
+{
+ Variable y= Variable (2);
+ Variable x= Variable (1);
+
+ nmod_poly_t f;
+ mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
+ nmod_poly_init_preinv (f, getCharacteristic(), ninv);
+ nmod_poly_set (f, F);
+
+ nmod_poly_t buf;
+ CanonicalForm result= 0;
+ int i= 0;
+ int degf= nmod_poly_degree(f);
+ int k= 0;
+ int degfSubK, repLength, j;
+ while (degf >= k)
+ {
+ degfSubK= degf - k;
+ if (degfSubK >= d)
+ repLength= d;
+ else
+ repLength= degfSubK + 1;
+
+ nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
+ for (j= 0; j < repLength; j++)
+ nmod_poly_set_coeff_ui (buf, j, nmod_poly_get_coeff_ui (f, j + k));
+ _nmod_poly_normalise (buf);
+
+ result += convertnmod_poly_t2FacCF (buf, x)*power (y, i);
+ i++;
+ k= d*i;
+ nmod_poly_clear (buf);
+ }
+ nmod_poly_clear (f);
+
+ return result;
+}
+
+CanonicalForm
+mulMod2FLINTFpReci (const CanonicalForm& F, const CanonicalForm& G, const
+ CanonicalForm& M)
+{
+ int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
+ d1 /= 2;
+ d1 += 1;
+
+ nmod_poly_t F1, F2;
+ mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
+ nmod_poly_init_preinv (F1, getCharacteristic(), ninv);
+ nmod_poly_init_preinv (F2, getCharacteristic(), ninv);
+ kronSubRecipro (F1, F2, F, d1);
+
+ nmod_poly_t G1, G2;
+ nmod_poly_init_preinv (G1, getCharacteristic(), ninv);
+ nmod_poly_init_preinv (G2, getCharacteristic(), ninv);
+ kronSubRecipro (G1, G2, G, d1);
+
+ int k= d1*degree (M);
+ nmod_poly_mullow (F1, F1, G1, (long) k);
+
+ int degtailF= degree (tailcoeff (F), 1);;
+ int degtailG= degree (tailcoeff (G), 1);
+ int taildegF= taildegree (F);
+ int taildegG= taildegree (G);
+
+ int b= nmod_poly_degree (F2) + nmod_poly_degree (G2) - k - degtailF - degtailG
+ + d1*(2+taildegF + taildegG);
+ nmod_poly_mulhigh (F2, F2, G2, b);
+ nmod_poly_shift_right (F2, F2, b);
+ int d2= tmax (nmod_poly_degree (F2)/d1, nmod_poly_degree (F1)/d1);
+
+
+ CanonicalForm result= reverseSubst (F1, F2, d1, d2);
+
+ nmod_poly_clear (F1);
+ nmod_poly_clear (F2);
+ nmod_poly_clear (G1);
+ nmod_poly_clear (G2);
+ return result;
+}
+
+CanonicalForm
+mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
+ CanonicalForm& M)
+{
+ CanonicalForm A= F;
+ CanonicalForm B= G;
+
+ int degAx= degree (A, 1);
+ int degAy= degree (A, 2);
+ int degBx= degree (B, 1);
+ int degBy= degree (B, 2);
+ int d1= degAx + 1 + degBx;
+ int d2= tmax (degAy, degBy);
+
+ if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
+ return mulMod2FLINTFpReci (A, B, M);
+
+ nmod_poly_t FLINTA, FLINTB;
+ mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
+ nmod_poly_init_preinv (FLINTA, getCharacteristic(), ninv);
+ nmod_poly_init_preinv (FLINTB, getCharacteristic(), ninv);
+ kronSubFp (FLINTA, A, d1);
+ kronSubFp (FLINTB, B, d1);
+
+ int k= d1*degree (M);
+ nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
+
+ A= reverseSubstFp (FLINTA, d1);
+
+ nmod_poly_clear (FLINTA);
+ nmod_poly_clear (FLINTB);
+ return A;
+}
+
+CanonicalForm
+mulMod2FLINTQReci (const CanonicalForm& F, const CanonicalForm& G, const
+ CanonicalForm& M)
+{
+ int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
+ d1 /= 2;
+ d1 += 1;
+
+ fmpz_poly_t F1, F2;
+ fmpz_poly_init (F1);
+ fmpz_poly_init (F2);
+ kronSubRecipro (F1, F2, F, d1);
+
+ fmpz_poly_t G1, G2;
+ fmpz_poly_init (G1);
+ fmpz_poly_init (G2);
+ kronSubRecipro (G1, G2, G, d1);
+
+ int k= d1*degree (M);
+ fmpz_poly_mullow (F1, F1, G1, (long) k);
+
+ int degtailF= degree (tailcoeff (F), 1);;
+ int degtailG= degree (tailcoeff (G), 1);
+ int taildegF= taildegree (F);
+ int taildegG= taildegree (G);
+
+ int b= fmpz_poly_degree (F2) + fmpz_poly_degree (G2) - k - degtailF - degtailG
+ + d1*(2+taildegF + taildegG);
+ fmpz_poly_mulhigh_n (F2, F2, G2, b);
+ fmpz_poly_shift_right (F2, F2, b);
+ int d2= tmax (fmpz_poly_degree (F2)/d1, fmpz_poly_degree (F1)/d1);
+
+ CanonicalForm result= reverseSubst (F1, F2, d1, d2);
+
+ fmpz_poly_clear (F1);
+ fmpz_poly_clear (F2);
+ fmpz_poly_clear (G1);
+ fmpz_poly_clear (G2);
+ return result;
+}
+
+CanonicalForm
+mulMod2FLINTQ (const CanonicalForm& F, const CanonicalForm& G, const
+ CanonicalForm& M)
+{
+ CanonicalForm A= F;
+ CanonicalForm B= G;
+
+ int degAx= degree (A, 1);
+ int degBx= degree (B, 1);
+ int d1= degAx + 1 + degBx;
+
+ CanonicalForm f= bCommonDen (F);
+ CanonicalForm g= bCommonDen (G);
+ A *= f;
+ B *= g;
+
+ fmpz_poly_t FLINTA, FLINTB;
+ fmpz_poly_init (FLINTA);
+ fmpz_poly_init (FLINTB);
+ kronSub (FLINTA, A, d1);
+ kronSub (FLINTB, B, d1);
+ int k= d1*degree (M);
+
+ fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
+ A= reverseSubstQ (FLINTA, d1);
+ fmpz_poly_clear (FLINTA);
+ fmpz_poly_clear (FLINTB);
+ return A/(f*g);
+}
+
+#endif
+
+CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
+ const CanonicalForm& M)
+{
+ if (A.isZero() || B.isZero())
+ return 0;
+
+ ASSERT (M.isUnivariate(), "M must be univariate");
+
+ CanonicalForm F= mod (A, M);
+ CanonicalForm G= mod (B, M);
+ if (F.inCoeffDomain() || G.inCoeffDomain())
+ return F*G;
+ Variable y= M.mvar();
+ int degF= degree (F, y);
+ int degG= degree (G, y);
+
+ if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) &&
+ (F.level() == G.level()))
+ {
+ CanonicalForm result= mulNTL (F, G);
+ return mod (result, M);
+ }
+ else if (degF <= 1 && degG <= 1)
+ {
+ CanonicalForm result= F*G;
+ return mod (result, M);
+ }
+
+ int sizeF= size (F);
+ int sizeG= size (G);
+
+ int fallBackToNaive= 50;
+ if (sizeF < fallBackToNaive || sizeG < fallBackToNaive)
+ return mod (F*G, M);
+
+#ifdef HAVE_FLINT
+ Variable alpha;
+ if (getCharacteristic() == 0 && !hasFirstAlgVar (F, alpha) && ! hasFirstAlgVar (G, alpha))
+ return mulMod2FLINTQ (F, G, M);
+#endif
+
+ if (getCharacteristic() > 0 && CFFactory::gettype() != GaloisFieldDomain &&
+ (((degF-degG) < 50 && degF > degG) || ((degG-degF) < 50 && degF <= degG)))
+ return mulMod2NTLFq (F, G, M);
+
+ int m= (int) ceil (degree (M)/2.0);
+ if (degF >= m || degG >= m)
+ {
+ CanonicalForm MLo= power (y, m);
+ CanonicalForm MHi= power (y, degree (M) - m);
+ CanonicalForm F0= mod (F, MLo);
+ CanonicalForm F1= div (F, MLo);
+ CanonicalForm G0= mod (G, MLo);
+ CanonicalForm G1= div (G, MLo);
+ CanonicalForm F0G1= mulMod2 (F0, G1, MHi);
+ CanonicalForm F1G0= mulMod2 (F1, G0, MHi);
+ CanonicalForm F0G0= mulMod2 (F0, G0, M);
+ return F0G0 + MLo*(F0G1 + F1G0);
+ }
+ else
+ {
+ m= (int) ceil (tmax (degF, degG)/2.0);
+ CanonicalForm yToM= power (y, m);
+ CanonicalForm F0= mod (F, yToM);
+ CanonicalForm F1= div (F, yToM);
+ CanonicalForm G0= mod (G, yToM);
+ CanonicalForm G1= div (G, yToM);
+ CanonicalForm H00= mulMod2 (F0, G0, M);
+ CanonicalForm H11= mulMod2 (F1, G1, M);
+ CanonicalForm H01= mulMod2 (F0 + F1, G0 + G1, M);
+ return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
+ }
+ DEBOUTLN (cerr, "fatal end in mulMod2");
+}
+
+CanonicalForm mod (const CanonicalForm& F, const CFList& M)
+{
+ CanonicalForm A= F;
+ for (CFListIterator i= M; i.hasItem(); i++)
+ A= mod (A, i.getItem());
+ return A;
+}
+
+CanonicalForm mulMod (const CanonicalForm& A, const CanonicalForm& B,
+ const CFList& MOD)
+{
+ if (A.isZero() || B.isZero())
+ return 0;
+
+ if (MOD.length() == 1)
+ return mulMod2 (A, B, MOD.getLast());
+
+ CanonicalForm M= MOD.getLast();
+ CanonicalForm F= mod (A, M);
+ CanonicalForm G= mod (B, M);
+ if (F.inCoeffDomain() || G.inCoeffDomain())
+ return F*G;
+ Variable y= M.mvar();
+ int degF= degree (F, y);
+ int degG= degree (G, y);
+
+ if ((degF <= 1 && F.level() <= M.level()) &&
+ (degG <= 1 && G.level() <= M.level()))
+ {
+ CFList buf= MOD;
+ buf.removeLast();
+ if (degF == 1 && degG == 1)
+ {
+ CanonicalForm F0= mod (F, y);
+ CanonicalForm F1= div (F, y);
+ CanonicalForm G0= mod (G, y);
+ CanonicalForm G1= div (G, y);
+ if (degree (M) > 2)
+ {
+ CanonicalForm H00= mulMod (F0, G0, buf);
+ CanonicalForm H11= mulMod (F1, G1, buf);
+ CanonicalForm H01= mulMod (F0 + F1, G0 + G1, buf);
+ return H11*y*y + (H01 - H00 - H11)*y + H00;
+ }
+ else //here degree (M) == 2
+ {
+ buf.append (y);
+ CanonicalForm F0G1= mulMod (F0, G1, buf);
+ CanonicalForm F1G0= mulMod (F1, G0, buf);
+ CanonicalForm F0G0= mulMod (F0, G0, MOD);
+ CanonicalForm result= F0G0 + y*(F0G1 + F1G0);
+ return result;
+ }
+ }
+ else if (degF == 1 && degG == 0)
+ return mulMod (div (F, y), G, buf)*y + mulMod (mod (F, y), G, buf);
+ else if (degF == 0 && degG == 1)
+ return mulMod (div (G, y), F, buf)*y + mulMod (mod (G, y), F, buf);
+ else
+ return mulMod (F, G, buf);
+ }
+ int m= (int) ceil (degree (M)/2.0);
+ if (degF >= m || degG >= m)
+ {
+ CanonicalForm MLo= power (y, m);
+ CanonicalForm MHi= power (y, degree (M) - m);
+ CanonicalForm F0= mod (F, MLo);
+ CanonicalForm F1= div (F, MLo);
+ CanonicalForm G0= mod (G, MLo);
+ CanonicalForm G1= div (G, MLo);
+ CFList buf= MOD;
+ buf.removeLast();
+ buf.append (MHi);
+ CanonicalForm F0G1= mulMod (F0, G1, buf);
+ CanonicalForm F1G0= mulMod (F1, G0, buf);
+ CanonicalForm F0G0= mulMod (F0, G0, MOD);
+ return F0G0 + MLo*(F0G1 + F1G0);
+ }
+ else
+ {
+ m= (int) ceil (tmax (degF, degG)/2.0);
+ CanonicalForm yToM= power (y, m);
+ CanonicalForm F0= mod (F, yToM);
+ CanonicalForm F1= div (F, yToM);
+ CanonicalForm G0= mod (G, yToM);
+ CanonicalForm G1= div (G, yToM);
+ CanonicalForm H00= mulMod (F0, G0, MOD);
+ CanonicalForm H11= mulMod (F1, G1, MOD);
+ CanonicalForm H01= mulMod (F0 + F1, G0 + G1, MOD);
+ return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
+ }
+ DEBOUTLN (cerr, "fatal end in mulMod");
+}
+
+CanonicalForm prodMod (const CFList& L, const CanonicalForm& M)
+{
+ if (L.isEmpty())
+ return 1;
+ int l= L.length();
+ if (l == 1)
+ return mod (L.getFirst(), M);
+ else if (l == 2) {
+ CanonicalForm result= mulMod2 (L.getFirst(), L.getLast(), M);
+ return result;
+ }
+ else
+ {
+ l /= 2;
+ CFList tmp1, tmp2;
+ CFListIterator i= L;
+ CanonicalForm buf1, buf2;
+ for (int j= 1; j <= l; j++, i++)
+ tmp1.append (i.getItem());
+ tmp2= Difference (L, tmp1);
+ buf1= prodMod (tmp1, M);
+ buf2= prodMod (tmp2, M);
+ CanonicalForm result= mulMod2 (buf1, buf2, M);
+ return result;
+ }
+}
+
+CanonicalForm prodMod (const CFList& L, const CFList& M)
+{
+ if (L.isEmpty())
+ return 1;
+ else if (L.length() == 1)
+ return L.getFirst();
+ else if (L.length() == 2)
+ return mulMod (L.getFirst(), L.getLast(), M);
+ else
+ {
+ int l= L.length()/2;
+ CFListIterator i= L;
+ CFList tmp1, tmp2;
+ CanonicalForm buf1, buf2;
+ for (int j= 1; j <= l; j++, i++)
+ tmp1.append (i.getItem());
+ tmp2= Difference (L, tmp1);
+ buf1= prodMod (tmp1, M);
+ buf2= prodMod (tmp2, M);
+ return mulMod (buf1, buf2, M);
+ }
+}
+
+CanonicalForm reverse (const CanonicalForm& F, int d)
+{
+ if (d == 0)
+ return F;
+ CanonicalForm A= F;
+ Variable y= Variable (2);
+ Variable x= Variable (1);
+ if (degree (A, x) > 0)
+ {
+ A= swapvar (A, x, y);
+ CanonicalForm result= 0;
+ CFIterator i= A;
+ while (d - i.exp() < 0)
+ i++;
+
+ for (; i.hasTerms() && (d - i.exp() >= 0); i++)
+ result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
+ return result;
+ }
+ else
+ return A*power (x, d);
+}
+
+CanonicalForm
+newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
+{
+ int l= ilog2(n);
+
+ CanonicalForm g= mod (F, M)[0] [0];
+
+ ASSERT (!g.isZero(), "expected a unit");
+
+ Variable alpha;
+
+ if (!g.isOne())
+ g = 1/g;
+ Variable x= Variable (1);
+ CanonicalForm result;
+ int exp= 0;
+ if (n & 1)
+ {
+ result= g;
+ exp= 1;
+ }
+ CanonicalForm h;
+
+ for (int i= 1; i <= l; i++)
+ {
+ h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
+ h= mod (h, power (x, (1 << i)) - 1);
+ h= div (h, power (x, (1 << (i - 1))));
+ h= mod (h, M);
+ g -= power (x, (1 << (i - 1)))*
+ mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
+
+ if (n & (1 << i))
+ {
+ if (exp)
+ {
+ h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
+ h= mod (h, power (x, exp + (1 << i)) - 1);
+ h= div (h, power (x, exp));
+ h= mod (h, M);
+ result -= power(x, exp)*mod (mulMod2 (g, h, M),
+ power (x, (1 << i)));
+ exp += (1 << i);
+ }
+ else
+ {
+ exp= (1 << i);
+ result= g;
+ }
+ }
+ }
+
+ return result;
+}
+
+CanonicalForm
+newtonDiv (const CanonicalForm& F, const CanonicalForm& G, const CanonicalForm&
+ M)
+{
+ ASSERT (getCharacteristic() > 0, "positive characteristic expected");
+ ASSERT (CFFactory::gettype() != GaloisFieldDomain, "no GF expected");
+
+ CanonicalForm A= mod (F, M);
+ CanonicalForm B= mod (G, M);
+
+ Variable x= Variable (1);
+ int degA= degree (A, x);
+ int degB= degree (B, x);
+ int m= degA - degB;
+ if (m < 0)
+ return 0;
+
+ Variable v;
+ CanonicalForm Q;
+ if (degB < 1 || CFFactory::gettype() == GaloisFieldDomain)
+ {
+ CanonicalForm R;
+ divrem2 (A, B, Q, R, M);
+ }
+ else
+ {
+ if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
+ {
+ CanonicalForm R= reverse (A, degA);
+ CanonicalForm revB= reverse (B, degB);
+ revB= newtonInverse (revB, m + 1, M);
+ Q= mulMod2 (R, revB, M);
+ Q= mod (Q, power (x, m + 1));
+ Q= reverse (Q, m);
+ }
+ else
+ {
+ zz_pX mipo= convertFacCF2NTLzzpX (M);
+ Variable y= Variable (2);
+ zz_pEX NTLA, NTLB;
+ NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
+ NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
+ div (NTLA, NTLA, NTLB);
+ Q= convertNTLzz_pEX2CF (NTLA, x, y);
+ }
+ }
+
+ return Q;
+}
+
+void
+newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+ CanonicalForm& R, const CanonicalForm& M)
+{
+ CanonicalForm A= mod (F, M);
+ CanonicalForm B= mod (G, M);
+ Variable x= Variable (1);
+ int degA= degree (A, x);
+ int degB= degree (B, x);
+ int m= degA - degB;
+
+ if (m < 0)
+ {
+ R= A;
+ Q= 0;
+ return;
+ }
+
+ Variable v;
+ if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
+ {
+ divrem2 (A, B, Q, R, M);
+ }
+ else
+ {
+ if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
+ {
+ R= reverse (A, degA);
+
+ CanonicalForm revB= reverse (B, degB);
+ revB= newtonInverse (revB, m + 1, M);
+ Q= mulMod2 (R, revB, M);
+
+ Q= mod (Q, power (x, m + 1));
+ Q= reverse (Q, m);
+
+ R= A - mulMod2 (Q, B, M);
+ }
+ else
+ {
+ zz_pX mipo= convertFacCF2NTLzzpX (M);
+ Variable y= Variable (2);
+ zz_pEX NTLA, NTLB;
+ NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
+ NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
+ zz_pEX NTLQ, NTLR;
+ DivRem (NTLQ, NTLR, NTLA, NTLB);
+ Q= convertNTLzz_pEX2CF (NTLQ, x, y);
+ R= convertNTLzz_pEX2CF (NTLR, x, y);
+ }
+ }
+}
+
+static inline
+CFList split (const CanonicalForm& F, const int m, const Variable& x)
+{
+ CanonicalForm A= F;
+ CanonicalForm buf= 0;
+ bool swap= false;
+ if (degree (A, x) <= 0)
+ return CFList(A);
+ else if (x.level() != A.level())
+ {
+ swap= true;
+ A= swapvar (A, x, A.mvar());
+ }
+
+ int j= (int) floor ((double) degree (A)/ m);
+ CFList result;
+ CFIterator i= A;
+ for (; j >= 0; j--)
+ {
+ while (i.hasTerms() && i.exp() - j*m >= 0)
+ {
+ if (swap)
+ buf += i.coeff()*power (A.mvar(), i.exp() - j*m);
+ else
+ buf += i.coeff()*power (x, i.exp() - j*m);
+ i++;
+ }
+ if (swap)
+ result.append (swapvar (buf, x, F.mvar()));
+ else
+ result.append (buf);
+ buf= 0;
+ }
+ return result;
+}
+
+static inline
+void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+ CanonicalForm& R, const CFList& M);
+
+static inline
+void divrem21 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+ CanonicalForm& R, const CFList& M)
+{
+ CanonicalForm A= mod (F, M);
+ CanonicalForm B= mod (G, M);
+ Variable x= Variable (1);
+ int degB= degree (B, x);
+ int degA= degree (A, x);
+ if (degA < degB)
+ {
+ Q= 0;
+ R= A;
+ return;
+ }
+ ASSERT (2*degB > degA, "expected degree (F, 1) < 2*degree (G, 1)");
+ if (degB < 1)
+ {
+ divrem (A, B, Q, R);
+ Q= mod (Q, M);
+ R= mod (R, M);
+ return;
+ }
+
+ int m= (int) ceil ((double) (degB + 1)/2.0) + 1;
+ CFList splitA= split (A, m, x);
+ if (splitA.length() == 3)
+ splitA.insert (0);
+ if (splitA.length() == 2)
+ {
+ splitA.insert (0);
+ splitA.insert (0);
+ }
+ if (splitA.length() == 1)
+ {
+ splitA.insert (0);
+ splitA.insert (0);
+ splitA.insert (0);
+ }
+
+ CanonicalForm xToM= power (x, m);
+
+ CFListIterator i= splitA;
+ CanonicalForm H= i.getItem();
+ i++;
+ H *= xToM;
+ H += i.getItem();
+ i++;
+ H *= xToM;
+ H += i.getItem();
+ i++;
+
+ divrem32 (H, B, Q, R, M);
+
+ CFList splitR= split (R, m, x);
+ if (splitR.length() == 1)
+ splitR.insert (0);
+
+ H= splitR.getFirst();
+ H *= xToM;
+ H += splitR.getLast();
+ H *= xToM;
+ H += i.getItem();
+
+ CanonicalForm bufQ;
+ divrem32 (H, B, bufQ, R, M);
+
+ Q *= xToM;
+ Q += bufQ;
+ return;
+}
+
+static inline
+void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+ CanonicalForm& R, const CFList& M)
+{
+ CanonicalForm A= mod (F, M);
+ CanonicalForm B= mod (G, M);
+ Variable x= Variable (1);
+ int degB= degree (B, x);
+ int degA= degree (A, x);
+ if (degA < degB)
+ {
+ Q= 0;
+ R= A;
+ return;
+ }
+ ASSERT (3*(degB/2) > degA, "expected degree (F, 1) < 3*(degree (G, 1)/2)");
+ if (degB < 1)
+ {
+ divrem (A, B, Q, R);
+ Q= mod (Q, M);
+ R= mod (R, M);
+ return;
+ }
+ int m= (int) ceil ((double) (degB + 1)/ 2.0);
+
+ CFList splitA= split (A, m, x);
+ CFList splitB= split (B, m, x);
+
+ if (splitA.length() == 2)
+ {
+ splitA.insert (0);
+ }
+ if (splitA.length() == 1)
+ {
+ splitA.insert (0);
+ splitA.insert (0);
+ }
+ CanonicalForm xToM= power (x, m);
+
+ CanonicalForm H;
+ CFListIterator i= splitA;
+ i++;
+
+ if (degree (splitA.getFirst(), x) < degree (splitB.getFirst(), x))
+ {
+ H= splitA.getFirst()*xToM + i.getItem();
+ divrem21 (H, splitB.getFirst(), Q, R, M);
+ }
+ else
+ {
+ R= splitA.getFirst()*xToM + i.getItem() + splitB.getFirst() -
+ splitB.getFirst()*xToM;
+ Q= xToM - 1;
+ }
+
+ H= mulMod (Q, splitB.getLast(), M);
+
+ R= R*xToM + splitA.getLast() - H;
+
+ while (degree (R, x) >= degB)
+ {
+ xToM= power (x, degree (R, x) - degB);
+ Q += LC (R, x)*xToM;
+ R -= mulMod (LC (R, x), B, M)*xToM;
+ Q= mod (Q, M);
+ R= mod (R, M);
+ }
+
+ return;
+}
+
+void divrem2 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+ CanonicalForm& R, const CanonicalForm& M)
+{
+ CanonicalForm A= mod (F, M);
+ CanonicalForm B= mod (G, M);
+
+ if (B.inCoeffDomain())
+ {
+ divrem (A, B, Q, R);
+ return;
+ }
+ if (A.inCoeffDomain() && !B.inCoeffDomain())
+ {
+ Q= 0;
+ R= A;
+ return;
+ }
+
+ if (B.level() < A.level())
+ {
+ divrem (A, B, Q, R);
+ return;
+ }
+ if (A.level() > B.level())
+ {
+ R= A;
+ Q= 0;
+ return;
+ }
+ if (B.level() == 1 && B.isUnivariate())
+ {
+ divrem (A, B, Q, R);
+ return;
+ }
+ if (!(B.level() == 1 && B.isUnivariate()) && (A.level() == 1 && A.isUnivariate()))
+ {
+ Q= 0;
+ R= A;
+ return;
+ }
+
+ Variable x= Variable (1);
+ int degB= degree (B, x);
+ if (degB > degree (A, x))
+ {
+ Q= 0;
+ R= A;
+ return;
+ }
+
+ CFList splitA= split (A, degB, x);
+
+ CanonicalForm xToDegB= power (x, degB);
+ CanonicalForm H, bufQ;
+ Q= 0;
+ CFListIterator i= splitA;
+ H= i.getItem()*xToDegB;
+ i++;
+ H += i.getItem();
+ CFList buf;
+ while (i.hasItem())
+ {
+ buf= CFList (M);
+ divrem21 (H, B, bufQ, R, buf);
+ i++;
+ if (i.hasItem())
+ H= R*xToDegB + i.getItem();
+ Q *= xToDegB;
+ Q += bufQ;
+ }
+ return;
+}
+
+void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
+ CanonicalForm& R, const CFList& MOD)
+{
+ CanonicalForm A= mod (F, MOD);
+ CanonicalForm B= mod (G, MOD);
+ Variable x= Variable (1);
+ int degB= degree (B, x);
+ if (degB > degree (A, x))
+ {
+ Q= 0;
+ R= A;
+ return;
+ }
+
+ if (degB <= 0)
+ {
+ divrem (A, B, Q, R);
+ Q= mod (Q, MOD);
+ R= mod (R, MOD);
+ return;
+ }
+ CFList splitA= split (A, degB, x);
+
+ CanonicalForm xToDegB= power (x, degB);
+ CanonicalForm H, bufQ;
+ Q= 0;
+ CFListIterator i= splitA;
+ H= i.getItem()*xToDegB;
+ i++;
+ H += i.getItem();
+ while (i.hasItem())
+ {
+ divrem21 (H, B, bufQ, R, MOD);
+ i++;
+ if (i.hasItem())
+ H= R*xToDegB + i.getItem();
+ Q *= xToDegB;
+ Q += bufQ;
+ }
+ return;
+}
+
+#endif
diff --git a/factory/facMul.h b/factory/facMul.h
new file mode 100644
index 0000000..76337e0
--- /dev/null
+++ b/factory/facMul.h
@@ -0,0 +1,147 @@
+/*****************************************************************************\
+ * Computer Algebra System SINGULAR
+\*****************************************************************************/
+/** @file facHensel.h
+ *
+ * This file defines functions for fast multiplication and division with
+ * remainder
+ *
+ * @author Martin Lee
+ *
+ **/
+/*****************************************************************************/
+
+#ifndef FAC_MUL_H
+#define FAC_MUL_H
+
+#include "canonicalform.h"
+
+/// multiplication of univariate polys over a finite field using NTL, if we are
+/// in GF factory's default multiplication is used.
+///
+/// @return @a mulNTL returns F*G
+CanonicalForm
+mulNTL (const CanonicalForm& F, ///< [in] a univariate poly
+ const CanonicalForm& G ///< [in] a univariate poly
+ );
+
+/// mod of univariate polys over a finite field using NTL, if we are
+/// in GF factory's default mod is used.
+///
+/// @return @a modNTL returns F mod G
+CanonicalForm
+modNTL (const CanonicalForm& F, ///< [in] a univariate poly
+ const CanonicalForm& G ///< [in] a univariate poly
+ );
+
+/// division of univariate polys over a finite field using NTL, if we are
+/// in GF factory's default division is used.
+///
+/// @return @a divNTL returns F/G
+CanonicalForm
+divNTL (const CanonicalForm& F, ///< [in] a univariate poly
+ const CanonicalForm& G ///< [in] a univariate poly
+ );
+
+/// division with remainder of @a F by
+/// @a G wrt Variable (1) modulo @a M.
+///
+/// @return @a Q returns the dividend, @a R returns the remainder.
+/// @sa divrem()
+void divrem2 (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
+ const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
+ CanonicalForm& Q, ///< [in,out] dividend
+ CanonicalForm& R, ///< [in,out] remainder, degree (R, 1) <
+ ///< degree (G, 1)
+ const CanonicalForm& M ///< [in] power of Variable (2)
+ );
+
+/// division with remainder of @a F by
+/// @a G wrt Variable (1) modulo @a MOD.
+///
+/// @sa divrem2()
+void divrem (
+ const CanonicalForm& F, ///< [in] multivariate, compressed polynomial
+ const CanonicalForm& G, ///< [in] multivariate, compressed polynomial
+ CanonicalForm& Q, ///< [in,out] dividend
+ CanonicalForm& R, ///< [in,out] remainder, degree (R, 1) <
+ ///< degree (G, 1)
+ const CFList& MOD ///< [in] only contains powers of
+ ///< Variables of level higher than 1
+ );
+
+
+/// division with remainder of @a F by
+/// @a G wrt Variable (1) modulo @a M using Newton inversion
+///
+/// @return @a Q returns the dividend, @a R returns the remainder.
+/// @sa divrem2(), newtonDiv()
+void
+newtonDivrem (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
+ const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
+ ///< which is monic in Variable (1)
+ CanonicalForm& Q, ///< [in,out] dividend
+ CanonicalForm& R, ///< [in,out] remainder, degree (R, 1) <
+ ///< degree (G, 1)
+ const CanonicalForm& M ///< [in] power of Variable (2)
+ );
+
+/// division of @a F by
+/// @a G wrt Variable (1) modulo @a M using Newton inversion
+///
+/// @return @a newtonDiv returns the dividend
+/// @sa divrem2(), newtonDivrem()
+CanonicalForm
+newtonDiv (const CanonicalForm& F, ///< [in] bivariate, compressed polynomial
+ const CanonicalForm& G, ///< [in] bivariate, compressed polynomial
+ ///< which is monic in Variable (1)
+ const CanonicalForm& M ///< [in] power of Variable (2)
+ );
+
+/// Karatsuba style modular multiplication for bivariate polynomials.
+///
+/// @return @a mulMod2 returns @a A * @a B mod @a M.
+CanonicalForm
+mulMod2 (const CanonicalForm& A, ///< [in] bivariate, compressed polynomial
+ const CanonicalForm& B, ///< [in] bivariate, compressed polynomial
+ const CanonicalForm& M ///< [in] power of Variable (2)
+ );
+
+/// Karatsuba style modular multiplication for multivariate polynomials.
+///
+/// @return @a mulMod2 returns @a A * @a B mod @a MOD.
+CanonicalForm
+mulMod (const CanonicalForm& A, ///< [in] multivariate, compressed polynomial
+ const CanonicalForm& B, ///< [in] multivariate, compressed polynomial
+ const CFList& MOD ///< [in] only contains powers of
+ ///< Variables of level higher than 1
+ );
+
+/// reduce @a F modulo elements in @a M.
+///
+/// @return @a mod returns @a F modulo @a M
+CanonicalForm mod (const CanonicalForm& F, ///< [in] compressed polynomial
+ const CFList& M ///< [in] list containing only
+ ///< univariate polynomials
+ );
+
+/// product of all elements in @a L modulo @a M via divide-and-conquer.
+///
+/// @return @a prodMod returns product of all elements in @a L modulo @a M.
+CanonicalForm
+prodMod (const CFList& L, ///< [in] contains only bivariate, compressed
+ ///< polynomials
+ const CanonicalForm& M ///< [in] power of Variable (2)
+ );
+
+/// product of all elements in @a L modulo @a M via divide-and-conquer.
+///
+/// @return @a prodMod returns product of all elements in @a L modulo @a M.
+CanonicalForm
+prodMod (const CFList& L, ///< [in] contains multivariate, compressed
+ ///< polynomials
+ const CFList& M ///< [in] contains only powers of Variables
+ );
+#endif
+/* FAC_MUL_H */
+
--
an open source computer algebra system
More information about the debian-science-commits
mailing list