[SCM] an open source computer algebra system branch, cleanedupstream, updated. 6125e540ca6d66c307958938a9d53b245507c323
Bernhard R. Link
brlink at debian.org
Tue Apr 24 15:54:48 UTC 2012
The following commit has been merged in the cleanedupstream branch:
commit 04290e40560003ed5ac61a60f046642627c79e35
Author: Martin Lee <martinlee84 at web.de>
Date: Thu Feb 9 17:59:15 2012 +0100
chg: more reorganization
diff --git a/factory/cf_gcd_smallp.cc b/factory/cf_gcd_smallp.cc
index d4e6bc9..e79b205 100755
--- a/factory/cf_gcd_smallp.cc
+++ b/factory/cf_gcd_smallp.cc
@@ -22,6 +22,7 @@
#include "timing.h"
#include "canonicalform.h"
+#include "algext.h"
#include "cf_map.h"
#include "cf_util.h"
#include "templates/ftmpl_functions.h"
@@ -4026,8 +4027,8 @@ int Hensel_P (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
CFArray lcs= CFArray (2);
lcs [0]= shiftedLCsEval1.getFirst();
lcs [1]= shiftedLCsEval2.getFirst();
- henselLift122 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
- lcs, false);
+ nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
+ lcs, false);
for (CFListIterator i= factors; i.hasItem(); i++)
{
@@ -4043,9 +4044,9 @@ int Hensel_P (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
liftBounds[0]= liftBound;
for (int i= 1; i < U.level() - 1; i++)
liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
- factors= henselLift2 (UEval, factors, liftBounds, U.level() - 1, false,
- shiftedLCsEval1, shiftedLCsEval2, Pi, diophant,
- noOneToOne);
+ factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
+ false, shiftedLCsEval1, shiftedLCsEval2, Pi,
+ diophant, noOneToOne);
delete [] liftBounds;
if (noOneToOne)
return 0;
diff --git a/factory/facFactorize.cc b/factory/facFactorize.cc
index e75e731..a7b94c0 100644
--- a/factory/facFactorize.cc
+++ b/factory/facFactorize.cc
@@ -477,8 +477,8 @@ precomputeLeadingCoeff (const CanonicalForm& LCF, const CFList& LCFFactors,
CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
CFArray Pi;
CFList diophant;
- henselLift122 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
- leadingCoeffs, false);
+ nonMonicHenselLift12 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
+ leadingCoeffs, false);
if (sqrfPartF.level() > 2)
{
diff --git a/factory/facFqBivarUtil.cc b/factory/facFqBivarUtil.cc
index 9f48bd4..a596281 100644
--- a/factory/facFqBivarUtil.cc
+++ b/factory/facFqBivarUtil.cc
@@ -13,6 +13,7 @@
#include <config.h>
#include "cf_map.h"
+#include "algext.h"
#include "cf_map_ext.h"
#include "templates/ftmpl_functions.h"
#include "ExtensionInfo.h"
diff --git a/factory/facFqFactorize.cc b/factory/facFqFactorize.cc
index e6351d0..dacc4c9 100644
--- a/factory/facFqFactorize.cc
+++ b/factory/facFqFactorize.cc
@@ -1655,8 +1655,8 @@ precomputeLeadingCoeff (const CanonicalForm& LCF, const CFList& LCFFactors,
CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
CFArray Pi;
CFList diophant;
- henselLift122 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
- leadingCoeffs, false);
+ nonMonicHenselLift12 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
+ leadingCoeffs, false);
if (sqrfPartF.level() > 2)
{
diff --git a/factory/facHensel.cc b/factory/facHensel.cc
index bde2c1c..22fe683 100644
--- a/factory/facHensel.cc
+++ b/factory/facHensel.cc
@@ -20,9 +20,9 @@
#include "debug.h"
#include "timing.h"
+#include "algext.h"
#include "facHensel.h"
#include "facMul.h"
-#include "cf_util.h"
#include "fac_util.h"
#include "cf_algorithm.h"
#include "cf_primes.h"
@@ -120,7 +120,7 @@ void tryDiophantine (CFList& result, const CanonicalForm& F,
}
}
-static inline
+static
CFList mapinto (const CFList& L)
{
CFList result;
@@ -129,7 +129,7 @@ CFList mapinto (const CFList& L)
return result;
}
-static inline
+static
int mod (const CFList& L, const CanonicalForm& p)
{
for (CFListIterator i= L; i.hasItem(); i++)
@@ -141,7 +141,7 @@ int mod (const CFList& L, const CanonicalForm& p)
}
-static inline void
+static void
chineseRemainder (const CFList & x1, const CanonicalForm & q1,
const CFList & x2, const CanonicalForm & q2,
CFList & xnew, CanonicalForm & qnew)
@@ -157,7 +157,7 @@ chineseRemainder (const CFList & x1, const CanonicalForm & q1,
qnew= tmp2;
}
-static inline
+static
CFList Farey (const CFList& L, const CanonicalForm& q)
{
CFList result;
@@ -166,7 +166,7 @@ CFList Farey (const CFList& L, const CanonicalForm& q)
return result;
}
-static inline
+static
CFList replacevar (const CFList& L, const Variable& a, const Variable& b)
{
CFList result;
@@ -354,7 +354,6 @@ void sortList (CFList& list, const Variable& x)
}
}
-static inline
CFList diophantine (const CanonicalForm& F, const CFList& factors)
{
if (getCharacteristic() == 0)
@@ -664,9 +663,8 @@ henselLiftResume12 (const CanonicalForm& F, CFList& factors, int start, int
return;
}
-static inline
CFList
-biDiophantine (const CanonicalForm& F, const CFList& factors, const int d)
+biDiophantine (const CanonicalForm& F, const CFList& factors, int d)
{
Variable y= F.mvar();
CFList result;
@@ -765,10 +763,9 @@ biDiophantine (const CanonicalForm& F, const CFList& factors, const int d)
}
}
-static inline
CFList
multiRecDiophantine (const CanonicalForm& F, const CFList& factors,
- const CFList& recResult, const CFList& M, const int d)
+ const CFList& recResult, const CFList& M, int d)
{
Variable y= F.mvar();
CFList result;
@@ -854,7 +851,7 @@ multiRecDiophantine (const CanonicalForm& F, const CFList& factors,
return result;
}
-static inline
+static
void
henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
const CFList& diophant, CFMatrix& M, CFArray& Pi, int j,
@@ -1096,7 +1093,7 @@ henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
}
CFList
-henselLift23 (const CFList& eval, const CFList& factors, const int* l, CFList&
+henselLift23 (const CFList& eval, const CFList& factors, int* l, CFList&
diophant, CFArray& Pi, CFMatrix& M)
{
CFList buf= factors;
@@ -1163,8 +1160,7 @@ henselLiftResume (const CanonicalForm& F, CFList& factors, int start, int end,
CFList
henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
- diophant, CFArray& Pi, CFMatrix& M, const int lOld, const int
- lNew)
+ diophant, CFArray& Pi, CFMatrix& M, int lOld, int lNew)
{
diophant= multiRecDiophantine (F.getFirst(), factors, diophant, MOD, lOld);
int k= 0;
@@ -1204,8 +1200,8 @@ henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
}
CFList
-henselLift (const CFList& eval, const CFList& factors, const int* l, const int
- lLength, bool sort)
+henselLift (const CFList& eval, const CFList& factors, int* l, int lLength,
+ bool sort)
{
CFList diophant;
CFList buf= factors;
@@ -1238,10 +1234,12 @@ henselLift (const CFList& eval, const CFList& factors, const int* l, const int
return result;
}
+// nonmonic
+
void
-henselStep122 (const CanonicalForm& F, const CFList& factors,
- CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
- CFArray& Pi, int j, const CFArray& LCs)
+nonMonicHenselStep12 (const CanonicalForm& F, const CFList& factors,
+ CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
+ CFArray& Pi, int j, const CFArray& LCs)
{
Variable x= F.mvar();
CanonicalForm xToJ= power (x, j);
@@ -1450,8 +1448,9 @@ henselStep122 (const CanonicalForm& F, const CFList& factors,
}
void
-henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
- CFList& diophant, CFMatrix& M, const CFArray& LCs, bool sort)
+nonMonicHenselLift12 (const CanonicalForm& F, CFList& factors, int l,
+ CFArray& Pi, CFList& diophant, CFMatrix& M,
+ const CFArray& LCs, bool sort)
{
if (sort)
sortList (factors, Variable (1));
@@ -1517,7 +1516,7 @@ henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
}
for (i= 1; i < l; i++)
- henselStep122 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
+ nonMonicHenselStep12 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
factors= CFList();
for (i= 0; i < bufFactors.size(); i++)
@@ -1528,7 +1527,6 @@ henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
/// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$
/// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$
-static inline
CFList
diophantine (const CFList& recResult, const CFList& factors,
const CFList& products, const CFList& M, const CanonicalForm& E,
@@ -1607,11 +1605,11 @@ diophantine (const CFList& recResult, const CFList& factors,
return result;
}
-static inline
void
-henselStep2 (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
- const CFList& diophant, CFMatrix& M, CFArray& Pi,
- const CFList& products, int j, const CFList& MOD, bool& noOneToOne)
+nonMonicHenselStep (const CanonicalForm& F, const CFList& factors,
+ CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
+ CFArray& Pi, const CFList& products, int j,
+ const CFList& MOD, bool& noOneToOne)
{
CanonicalForm E;
CanonicalForm xToJ= power (F.mvar(), j);
@@ -1847,9 +1845,9 @@ CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c)
}
CFList
-henselLift232 (const CFList& eval, const CFList& factors, int* l, CFList&
- diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
- const CFList& LCs2, bool& bad)
+nonMonicHenselLift232(const CFList& eval, const CFList& factors, int* l, CFList&
+ diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
+ const CFList& LCs2, bool& bad)
{
CFList buf= factors;
int k= 0;
@@ -1897,7 +1895,8 @@ henselLift232 (const CFList& eval, const CFList& factors, int* l, CFList&
for (int d= 1; d < l[1]; d++)
{
- henselStep2 (j.getItem(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
+ nonMonicHenselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, products,
+ d, MOD, bad);
if (bad)
return CFList();
}
@@ -1909,9 +1908,10 @@ henselLift232 (const CFList& eval, const CFList& factors, int* l, CFList&
CFList
-henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
- diophant, CFArray& Pi, CFMatrix& M, const int lOld, int&
- lNew, const CFList& LCs1, const CFList& LCs2, bool& bad)
+nonMonicHenselLift2 (const CFList& F, const CFList& factors, const CFList& MOD,
+ CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
+ int& lNew, const CFList& LCs1, const CFList& LCs2, bool& bad
+ )
{
int k= 0;
CFArray bufFactors= CFArray (factors.length());
@@ -1958,7 +1958,8 @@ henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
for (int d= 1; d < lNew; d++)
{
- henselStep2 (F.getLast(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
+ nonMonicHenselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, products,
+ d, MOD, bad);
if (bad)
return CFList();
}
@@ -1970,9 +1971,9 @@ henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
}
CFList
-henselLift2 (const CFList& eval, const CFList& factors, int* l, const int
- lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
- const CFArray& Pi, const CFList& diophant, bool& bad)
+nonMonicHenselLift2 (const CFList& eval, const CFList& factors, int* l, int
+ lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
+ const CFArray& Pi, const CFList& diophant, bool& bad)
{
CFList bufDiophant= diophant;
CFList buf= factors;
@@ -1980,8 +1981,8 @@ henselLift2 (const CFList& eval, const CFList& factors, int* l, const int
sortList (buf, Variable (1));
CFArray bufPi= Pi;
CFMatrix M= CFMatrix (l[1], factors.length());
- CFList result= henselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2,
- bad);
+ CFList result=
+ nonMonicHenselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2, bad);
if (bad)
return CFList();
@@ -2009,8 +2010,8 @@ henselLift2 (const CFList& eval, const CFList& factors, int* l, const int
bufLCs1.append (jj.getItem());
bufLCs2.append (jjj.getItem());
M= CFMatrix (l[i], factors.length());
- result= henselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M, l[i - 1],
- l[i], bufLCs1, bufLCs2, bad);
+ result= nonMonicHenselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M,
+ l[i - 1], l[i], bufLCs1, bufLCs2, bad);
if (bad)
return CFList();
MOD.append (power (Variable (i + 2), l[i]));
@@ -2107,7 +2108,8 @@ nonMonicHenselLift23 (const CanonicalForm& F, const CFList& factors, const
MOD.insert (yToL);
for (int d= 1; d < liftBound; d++)
{
- henselStep2 (F, factors, bufFactors, diophant, M, Pi, products, d, MOD, bad);
+ nonMonicHenselStep (F, factors, bufFactors, diophant, M, Pi, products, d,
+ MOD, bad);
if (bad)
return CFList();
}
@@ -2120,7 +2122,7 @@ nonMonicHenselLift23 (const CanonicalForm& F, const CFList& factors, const
CFList
nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs,
- CFList& diophant, CFArray& Pi, CFMatrix& M, const int lOld,
+ CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
int& lNew, const CFList& MOD, bool& noOneToOne
)
{
@@ -2187,8 +2189,8 @@ nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs,
for (int d= 1; d < lNew; d++)
{
- henselStep2 (F.getLast(), factors, bufFactors, diophant, M, Pi, products, d,
- MOD, noOneToOne);
+ nonMonicHenselStep (F.getLast(), factors, bufFactors, diophant, M, Pi,
+ products, d, MOD, noOneToOne);
if (noOneToOne)
return CFList();
}
diff --git a/factory/facHensel.h b/factory/facHensel.h
index e3bf585..7cf857c 100644
--- a/factory/facHensel.h
+++ b/factory/facHensel.h
@@ -22,9 +22,6 @@
#include "timing.h"
#include "canonicalform.h"
-#include "cf_iter.h"
-#include "templates/ftmpl_functions.h"
-#include "algext.h"
/// sort a list of polynomials by their degree in @a x.
///
@@ -83,7 +80,7 @@ henselLift23 (const CFList& eval, ///< [in] contains compressed, bivariate
const CFList& factors, ///< [in] monic bivariate factors,
///< including leading coefficient
///< as first element.
- const int* l, ///< [in] l[0]: precision of bivariate
+ int* l, ///< [in] l[0]: precision of bivariate
///< lifting, l[1] as above
CFList& diophant, ///< [in,out] returns the result of
///< biDiophantine()
@@ -127,8 +124,8 @@ henselLift (const CFList& F, ///< [in] two compressed, multivariate
CFList& diophant, ///< [in,out] result of multiRecDiophantine()
CFArray& Pi, ///< [in,out] stores intermediate results
CFMatrix& M, ///< [in,out] stores intermediate results
- const int lOld, ///< [in] lifting precision of F.mvar()
- const int lNew ///< [in] lifting precision of G.mvar()
+ int lOld, ///< [in] lifting precision of F.mvar()
+ int lNew ///< [in] lifting precision of G.mvar()
);
/// Hensel lifting from bivariate to multivariate
@@ -144,16 +141,15 @@ henselLift (const CFList& eval, ///< [in] a list of polynomials the last
///< compressed bivariate poly.
const CFList& factors, ///< [in] bivariate factors, including
///< leading coefficient
- const int* l, ///< [in] lifting bounds
- const int lLength, ///< [in] length of l
+ int* l, ///< [in] lifting bounds
+ int lLength, ///< [in] length of l
bool sort= true ///< [in] sort factors by degree in
///< Variable(1)
);
-/// two factor Hensel lifting from univariate to bivariate, factors need not to
-/// be monic
+/// Hensel lifting from univariate to bivariate, factors need not to be monic
void
-henselLift122 (const CanonicalForm& F,///< [in] a bivariate poly
+nonMonicHenselLift12 (const CanonicalForm& F,///< [in] a bivariate poly
CFList& factors, ///< [in, out] a list of univariate polys
///< lifted factors
int l, ///< [in] lift bound
@@ -170,7 +166,7 @@ henselLift122 (const CanonicalForm& F,///< [in] a bivariate poly
///
/// @return @a henselLift122 returns a list of lifted factors
CFList
-henselLift2 (const CFList& eval, ///< [in] a list of polynomials the last
+nonMonicHenselLift2 (const CFList& eval,///< [in] a list of polynomials the last
///< element is a compressed multivariate
///< poly, last but one element equals the
///< last elements modulo its main variable
@@ -178,7 +174,7 @@ henselLift2 (const CFList& eval, ///< [in] a list of polynomials the last
///< compressed bivariate poly.
const CFList& factors,///< [in] bivariate factors
int* l, ///< [in] lift bounds
- const int lLength, ///< [in] length of l
+ int lLength, ///< [in] length of l
bool sort, ///< [in] if true factors are sorted by
///< their degree in Variable(1)
const CFList& LCs1, ///< [in] a list of evaluated LC of first
diff --git a/factory/facMul.cc b/factory/facMul.cc
index f022ac5..27e3b2f 100644
--- a/factory/facMul.cc
+++ b/factory/facMul.cc
@@ -16,6 +16,7 @@
#include "canonicalform.h"
#include "facMul.h"
#include "algext.h"
+#include "cf_util.h"
#include "templates/ftmpl_functions.h"
#ifdef HAVE_NTL
@@ -26,6 +27,8 @@
#include "FLINTconvert.h"
#endif
+// univariate polys
+
#ifdef HAVE_FLINT
void kronSub (fmpz_poly_t result, const CanonicalForm& A, int d)
{
@@ -47,8 +50,8 @@ void kronSub (fmpz_poly_t result, const CanonicalForm& A, int d)
CanonicalForm
-reverseSubst (const fmpz_poly_t F, int d, const Variable& alpha,
- const CanonicalForm& den)
+reverseSubstQa (const fmpz_poly_t F, int d, const Variable& alpha,
+ const CanonicalForm& den)
{
Variable x= Variable (1);
@@ -110,7 +113,7 @@ mulFLINTQa (const CanonicalForm& F, const CanonicalForm& G,
fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
denA *= denB;
- A= reverseSubst (FLINTA, d, alpha, denA);
+ A= reverseSubstQa (FLINTA, d, alpha, denA);
fmpz_poly_clear (FLINTA);
fmpz_poly_clear (FLINTB);
@@ -221,7 +224,7 @@ mulFLINTQaTrunc (const CanonicalForm& F, const CanonicalForm& G,
fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, k);
denA *= denB;
- A= reverseSubst (FLINTA, d, alpha, denA);
+ A= reverseSubstQa (FLINTA, d, alpha, denA);
fmpz_poly_clear (FLINTA);
fmpz_poly_clear (FLINTB);
return A;
@@ -569,6 +572,10 @@ divNTL (const CanonicalForm& F, const CanonicalForm& G)
return result;
}
+// end univariate polys
+//*************************
+// bivariate polys
+
#ifdef HAVE_FLINT
void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
{
@@ -591,34 +598,6 @@ void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
_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);
@@ -663,174 +642,149 @@ void kronSubQa (fmpq_poly_t result, const CanonicalForm& A, int d1, int d2)
fmpq_clear (coeff);
_fmpq_poly_normalise (result);
}
-#endif
-zz_pX kronSubFp (const CanonicalForm& A, int d)
+void
+kronSubReciproFp (nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm& A,
+ int d)
{
int degAy= degree (A);
- zz_pX result;
- result.rep.SetLength (d*(degAy + 1));
+ 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));
- zz_p *resultp;
- resultp= result.rep.elts();
- zz_pX buf;
- zz_p *bufp;
- int j, k, bufRepLength;
+ nmod_poly_t buf;
+ int k, kk, j, bufRepLength;
for (CFIterator i= A; i.hasTerms(); i++)
{
- if (i.coeff().inCoeffDomain())
- buf= convertFacCF2NTLzzpX (i.coeff());
- else
- buf= convertFacCF2NTLzzpX (i.coeff());
+ convertFacCF2nmod_poly_t (buf, i.coeff());
k= i.exp()*d;
- bufp= buf.rep.elts();
- bufRepLength= (int) buf.rep.length();
+ kk= (degAy - i.exp())*d;
+ bufRepLength= (int) nmod_poly_length (buf);
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));
+ 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()
+ )
+ );
}
- 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];
+ nmod_poly_clear (buf);
}
- result.normalize();
-
- return result;
+ _nmod_poly_normalise (subA1);
+ _nmod_poly_normalise (subA2);
}
void
-kronSubRecipro (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
- const Variable& alpha)
+kronSubReciproQ (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
+ int d)
{
int degAy= degree (A);
- subA1.rep.SetLength ((long) d*(degAy + 2));
- subA2.rep.SetLength ((long) d*(degAy + 2));
+ fmpz_poly_init2 (subA1, d*(degAy + 2));
+ fmpz_poly_init2 (subA2, 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;
+ fmpz_poly_t buf;
+ fmpz_t coeff1, coeff2;
+ int k, kk, j, 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);
+ convertFacCF2Fmpz_poly_t (buf, i.coeff());
k= i.exp()*d;
kk= (degAy - i.exp())*d;
- bufp= buf.rep.elts();
- bufRepLength= (int) buf.rep.length();
+ bufRepLength= (int) fmpz_poly_length (buf);
for (j= 0; j < bufRepLength; j++)
{
- subA1p [j + k] += bufp [j];
- subA2p [j + kk] += bufp [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);
}
- subA1.normalize();
- subA2.normalize();
+ fmpz_clear (coeff1);
+ fmpz_clear (coeff2);
+ _fmpz_poly_normalise (subA1);
+ _fmpz_poly_normalise (subA2);
}
-void
-kronSubRecipro (zz_pX& subA1, zz_pX& subA2,const CanonicalForm& A, int d)
+CanonicalForm reverseSubstQ (const fmpz_poly_t F, int d)
{
- int degAy= degree (A);
- subA1.rep.SetLength ((long) d*(degAy + 2));
- subA2.rep.SetLength ((long) d*(degAy + 2));
+ Variable y= Variable (2);
+ Variable x= Variable (1);
- 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;
+ fmpz_poly_t f;
+ fmpz_poly_init (f);
+ fmpz_poly_set (f, F);
- for (CFIterator i= A; i.hasTerms(); i++)
+ 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)
{
- buf= convertFacCF2NTLzzpX (i.coeff());
+ degfSubK= degf - k;
+ if (degfSubK >= d)
+ repLength= d;
+ else
+ repLength= degfSubK + 1;
- k= i.exp()*d;
- kk= (degAy - i.exp())*d;
- bufp= buf.rep.elts();
- bufRepLength= (int) buf.rep.length();
- for (j= 0; j < bufRepLength; j++)
+ fmpz_poly_init2 (buf, repLength);
+ fmpz_init (coeff);
+ for (j= 0; j < repLength; j++)
{
- subA1p [j + k] += bufp [j];
- subA2p [j + kk] += bufp [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);
}
- subA1.normalize();
- subA2.normalize();
+ fmpz_poly_clear (f);
+
+ return result;
}
CanonicalForm
-reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
- const Variable& alpha)
+reverseSubstReciproFp (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
{
Variable y= Variable (2);
Variable x= Variable (1);
- zz_pEX f= F;
- zz_pEX g= G;
- int degf= deg(f);
- int degg= deg(g);
+ 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);
- 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));
+ 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));
- zz_pE *gp= g.rep.elts();
- zz_pE *fp= f.rep.elts();
CanonicalForm result= 0;
int i= 0;
int lf= 0;
@@ -838,8 +792,6 @@ reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int 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)
@@ -848,14 +800,13 @@ reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
repLengthBuf1= 0;
else
repLengthBuf1= degfSubLf + 1;
- buf1.rep.SetLength((long) repLengthBuf1);
+ nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
- buf1p= buf1.rep.elts();
for (ind= 0; ind < repLengthBuf1; ind++)
- buf1p [ind]= fp [ind + lf];
- buf1.normalize();
+ nmod_poly_set_coeff_ui (buf1, ind, nmod_poly_get_coeff_ui (f, ind+lf));
+ _nmod_poly_normalise (buf1);
- repLengthBuf1= buf1.rep.length();
+ repLengthBuf1= nmod_poly_length (buf1);
if (deggSubLg >= d - 1)
repLengthBuf2= d - 1;
@@ -864,27 +815,23 @@ reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
else
repLengthBuf2= deggSubLg + 1;
- buf2.rep.SetLength ((long) repLengthBuf2);
- buf2p= buf2.rep.elts();
+ nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
for (ind= 0; ind < repLengthBuf2; ind++)
- buf2p [ind]= gp [ind + lg];
- buf2.normalize();
+ nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
- repLengthBuf2= buf2.rep.length();
+ _nmod_poly_normalise (buf2);
+ repLengthBuf2= nmod_poly_length (buf2);
- buf3.rep.SetLength((long) repLengthBuf2 + d);
- buf3p= buf3.rep.elts();
- buf2p= buf2.rep.elts();
- buf1p= buf1.rep.elts();
+ nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
for (ind= 0; ind < repLengthBuf1; ind++)
- buf3p [ind]= buf1p [ind];
+ nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
for (ind= repLengthBuf1; ind < d; ind++)
- buf3p [ind]= zzpEZero;
+ nmod_poly_set_coeff_ui (buf3, ind, 0);
for (ind= 0; ind < repLengthBuf2; ind++)
- buf3p [ind + d]= buf2p [ind];
- buf3.normalize();
+ nmod_poly_set_coeff_ui (buf3, ind+d, nmod_poly_get_coeff_ui (buf2, ind));
+ _nmod_poly_normalise (buf3);
- result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
+ result += convertnmod_poly_t2FacCF (buf3, x)*power (y, i);
i++;
@@ -894,55 +841,67 @@ reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
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];
+ 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;
-
- buf2p= buf2.rep.elts();
+ }
if (degfSubLf >= 0)
{
for (ind= 0; ind < repLengthBuf2; ind++)
- fp [ind + lf] -= buf2p [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 zz_pX& F, const zz_pX& G, int d, int k)
+reverseSubstReciproQ (const fmpz_poly_t F, const fmpz_poly_t 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);
+ 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);
- zz_pX buf1;
- zz_pX buf2;
- zz_pX buf3;
- zz_p *buf1p;
- zz_p *buf2p;
- zz_p *buf3p;
+ fmpz_poly_t buf1,buf2, buf3;
- if (f.rep.length() < (long) d*(k+1)) //zero padding
- f.rep.SetLength ((long)d*(k+1));
+ if (fmpz_poly_length (f) < (long) d*(k+1)) //zero padding
+ fmpz_poly_fit_length (f,(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;
@@ -950,7 +909,7 @@ reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
int degfSubLf= degf;
int deggSubLg= degg-lg;
int repLengthBuf2, repLengthBuf1, ind, tmp;
- zz_p zzpZero= zz_p();
+ fmpz_t tmp1, tmp2;
while (degf >= lf || lg >= 0)
{
if (degfSubLf >= d)
@@ -959,14 +918,17 @@ reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
repLengthBuf1= 0;
else
repLengthBuf1= degfSubLf + 1;
- buf1.rep.SetLength((long) repLengthBuf1);
- buf1p= buf1.rep.elts();
+ fmpz_poly_init2 (buf1, repLengthBuf1);
+
for (ind= 0; ind < repLengthBuf1; ind++)
- buf1p [ind]= fp [ind + lf];
- buf1.normalize();
+ {
+ fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
+ fmpz_poly_set_coeff_fmpz (buf1, ind, tmp1);
+ }
+ _fmpz_poly_normalise (buf1);
- repLengthBuf1= buf1.rep.length();
+ repLengthBuf1= fmpz_poly_length (buf1);
if (deggSubLg >= d - 1)
repLengthBuf2= d - 1;
@@ -975,29 +937,33 @@ reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
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();
+ fmpz_poly_init2 (buf2, repLengthBuf2);
- repLengthBuf2= buf2.rep.length();
+ 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);
- buf3.rep.SetLength((long) repLengthBuf2 + d);
- buf3p= buf3.rep.elts();
- buf2p= buf2.rep.elts();
- buf1p= buf1.rep.elts();
+ fmpz_poly_init2 (buf3, repLengthBuf2 + d);
for (ind= 0; ind < repLengthBuf1; ind++)
- buf3p [ind]= buf1p [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++)
- buf3p [ind]= zzpZero;
+ fmpz_poly_set_coeff_ui (buf3, ind, 0);
for (ind= 0; ind < repLengthBuf2; ind++)
- buf3p [ind + d]= buf2p [ind];
- buf3.normalize();
+ {
+ fmpz_poly_get_coeff_fmpz (tmp1, buf2, ind);
+ fmpz_poly_set_coeff_fmpz (buf3, ind + d, tmp1);
+ }
+ _fmpz_poly_normalise (buf3);
- result += convertNTLzzpX2CF (buf3, x)*power (y, i);
+ result += convertFmpz_poly_t2FacCF (buf3, x)*power (y, i);
i++;
@@ -1007,80 +973,63 @@ reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
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];
+ {
+ 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;
-
- buf2p= buf2.rep.elts();
+ }
if (degfSubLf >= 0)
{
for (ind= 0; ind < repLengthBuf2; ind++)
- fp [ind + lf] -= buf2p [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);
}
- 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;
- }
+ fmpz_poly_clear (f);
+ fmpz_poly_clear (g);
+ fmpz_clear (tmp1);
+ fmpz_clear (tmp2);
return result;
}
-CanonicalForm reverseSubstFp (const zz_pX& F, int d)
+CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
{
Variable y= Variable (2);
Variable x= Variable (1);
- zz_pX f= F;
- zz_p *fp= f.rep.elts();
+ nmod_poly_t f;
+ mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
+ nmod_poly_init_preinv (f, getCharacteristic(), ninv);
+ nmod_poly_set (f, F);
- zz_pX buf;
- zz_p *bufp;
+ nmod_poly_t buf;
CanonicalForm result= 0;
int i= 0;
- int degf= deg(f);
+ int degf= nmod_poly_degree(f);
int k= 0;
int degfSubK, repLength, j;
while (degf >= k)
@@ -1091,56 +1040,67 @@ CanonicalForm reverseSubstFp (const zz_pX& F, int d)
else
repLength= degfSubK + 1;
- buf.rep.SetLength ((long) repLength);
- bufp= buf.rep.elts();
+ nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
for (j= 0; j < repLength; j++)
- bufp [j]= fp [j + k];
- buf.normalize();
+ nmod_poly_set_coeff_ui (buf, j, nmod_poly_get_coeff_ui (f, j + k));
+ _nmod_poly_normalise (buf);
- result += convertNTLzzpX2CF (buf, x)*power (y, i);
+ result += convertnmod_poly_t2FacCF (buf, x)*power (y, i);
i++;
k= d*i;
+ nmod_poly_clear (buf);
}
+ nmod_poly_clear (f);
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)
+mulMod2FLINTFpReci (const CanonicalForm& F, const CanonicalForm& G, const
+ CanonicalForm& M)
{
- int d1= degree (F, 1) + degree (G, 1) + 1;
+ int d1= tmax (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);
+ 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);
+ kronSubReciproFp (F1, F2, F, d1);
+
+ nmod_poly_t G1, G2;
+ nmod_poly_init_preinv (G1, getCharacteristic(), ninv);
+ nmod_poly_init_preinv (G2, getCharacteristic(), ninv);
+ kronSubReciproFp (G1, G2, G, d1);
int k= d1*degree (M);
- MulTrunc (F1, F1, G1, (long) k);
+ nmod_poly_mullow (F1, F1, G1, (long) k);
- int degtailF= degree (tailcoeff (F), 1);
+ 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 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);
- int d2= tmax (deg (F2)/d1, deg (F1)/d1);
- return reverseSubst (F1, F2, d1, d2);
+
+ CanonicalForm result= reverseSubstReciproFp (F1, F2, d1, d2);
+
+ nmod_poly_clear (F1);
+ nmod_poly_clear (F2);
+ nmod_poly_clear (G1);
+ nmod_poly_clear (G2);
+ return result;
}
-//Kronecker substitution
CanonicalForm
-mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M)
+mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
+ CanonicalForm& M)
{
CanonicalForm A= F;
CanonicalForm B= G;
@@ -1153,244 +1113,264 @@ mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
int d2= tmax (degAy, degBy);
if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
- return mulMod2NTLFpReci (A, B, M);
+ return mulMod2FLINTFpReci (A, B, M);
- zz_pX NTLA= kronSubFp (A, d1);
- zz_pX NTLB= kronSubFp (B, d1);
+ 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);
- MulTrunc (NTLA, NTLA, NTLB, (long) k);
+ nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
- A= reverseSubstFp (NTLA, d1);
+ A= reverseSubstFp (FLINTA, d1);
+ nmod_poly_clear (FLINTA);
+ nmod_poly_clear (FLINTB);
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;
+mulMod2FLINTQReci (const CanonicalForm& F, const CanonicalForm& G, const
+ CanonicalForm& M)
+{
+ int d1= tmax (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);
+ fmpz_poly_t F1, F2;
+ fmpz_poly_init (F1);
+ fmpz_poly_init (F2);
+ kronSubReciproQ (F1, F2, F, d1);
+
+ fmpz_poly_t G1, G2;
+ fmpz_poly_init (G1);
+ fmpz_poly_init (G2);
+ kronSubReciproQ (G1, G2, G, d1);
int k= d1*degree (M);
- MulTrunc (F1, F1, G1, (long) k);
+ fmpz_poly_mullow (F1, F1, G1, (long) k);
- int degtailF= degree (tailcoeff (F), 1);
+ 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 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);
- int d2= tmax (deg (F2)/d1, deg (F1)/d1);
- return reverseSubst (F1, F2, d1, d2, alpha);
-}
+ CanonicalForm result= reverseSubstReciproQ (F1, F2, d1, d2);
-#ifdef HAVE_FLINT
-CanonicalForm
-mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M);
-#endif
+ fmpz_poly_clear (F1);
+ fmpz_poly_clear (F2);
+ fmpz_poly_clear (G1);
+ fmpz_poly_clear (G2);
+ return result;
+}
CanonicalForm
-mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M)
+mulMod2FLINTQ (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 degAx= degree (A, 1);
+ int degBx= degree (B, 1);
+ int d1= degAx + 1 + degBx;
- int k= d1*degree (M);
+ CanonicalForm f= bCommonDen (F);
+ CanonicalForm g= bCommonDen (G);
+ A *= f;
+ B *= g;
- MulTrunc (NTLA, NTLA, NTLB, (long) k);
+ 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);
- A= reverseSubst (NTLA, d1, alpha);
+ fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
+ A= reverseSubstQ (FLINTA, d1);
+ fmpz_poly_clear (FLINTA);
+ fmpz_poly_clear (FLINTB);
+ return A/(f*g);
+}
- 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)
+zz_pX kronSubFp (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));
+ zz_pX result;
+ result.rep.SetLength (d*(degAy + 1));
- nmod_poly_t buf;
+ zz_p *resultp;
+ resultp= result.rep.elts();
+ zz_pX buf;
+ zz_p *bufp;
+ int j, k, bufRepLength;
- int k, kk, j, bufRepLength;
for (CFIterator i= A; i.hasTerms(); i++)
{
- convertFacCF2nmod_poly_t (buf, i.coeff());
+ if (i.coeff().inCoeffDomain())
+ buf= convertFacCF2NTLzzpX (i.coeff());
+ else
+ buf= convertFacCF2NTLzzpX (i.coeff());
k= i.exp()*d;
- kk= (degAy - i.exp())*d;
- bufRepLength= (int) nmod_poly_length (buf);
+ 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 kronSubFq (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())
{
- 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()
- )
- );
+ buf2= convertFacCF2NTLzzpX (i.coeff());
+ buf1= to_zz_pEX (to_zz_pE (buf2));
}
- nmod_poly_clear (buf);
+ 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];
}
- _nmod_poly_normalise (subA1);
- _nmod_poly_normalise (subA2);
+ result.normalize();
+
+ return result;
}
void
-kronSubRecipro (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
- int d)
+kronSubReciproFq (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
+ const Variable& alpha)
{
int degAy= degree (A);
- fmpz_poly_init2 (subA1, d*(degAy + 2));
- fmpz_poly_init2 (subA2, d*(degAy + 2));
+ subA1.rep.SetLength ((long) d*(degAy + 2));
+ subA2.rep.SetLength ((long) d*(degAy + 2));
- fmpz_poly_t buf;
- fmpz_t coeff1, coeff2;
+ 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;
- int k, kk, j, bufRepLength;
for (CFIterator i= A; i.hasTerms(); i++)
{
- convertFacCF2Fmpz_poly_t (buf, i.coeff());
+ 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;
- bufRepLength= (int) fmpz_poly_length (buf);
+ bufp= buf.rep.elts();
+ bufRepLength= (int) buf.rep.length();
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);
+ subA1p [j + k] += bufp [j];
+ subA2p [j + kk] += bufp [j];
}
- fmpz_poly_clear (buf);
}
- fmpz_clear (coeff1);
- fmpz_clear (coeff2);
- _fmpz_poly_normalise (subA1);
- _fmpz_poly_normalise (subA2);
+ subA1.normalize();
+ subA2.normalize();
}
-CanonicalForm reverseSubstQ (const fmpz_poly_t F, int d)
+void
+kronSubReciproFp (zz_pX& subA1, zz_pX& subA2, const CanonicalForm& A, int d)
{
- Variable y= Variable (2);
- Variable x= Variable (1);
+ int degAy= degree (A);
+ subA1.rep.SetLength ((long) d*(degAy + 2));
+ subA2.rep.SetLength ((long) d*(degAy + 2));
- fmpz_poly_t f;
- fmpz_poly_init (f);
- fmpz_poly_set (f, F);
+ 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;
- 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)
+ for (CFIterator i= A; i.hasTerms(); i++)
{
- degfSubK= degf - k;
- if (degfSubK >= d)
- repLength= d;
- else
- repLength= degfSubK + 1;
+ buf= convertFacCF2NTLzzpX (i.coeff());
- fmpz_poly_init2 (buf, repLength);
- fmpz_init (coeff);
- for (j= 0; j < repLength; j++)
+ k= i.exp()*d;
+ kk= (degAy - i.exp())*d;
+ bufp= buf.rep.elts();
+ bufRepLength= (int) buf.rep.length();
+ for (j= 0; j < bufRepLength; j++)
{
- fmpz_poly_get_coeff_fmpz (coeff, f, j + k);
- fmpz_poly_set_coeff_fmpz (buf, j, coeff);
+ subA1p [j + k] += bufp [j];
+ subA2p [j + kk] += bufp [j];
}
- _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;
+ subA1.normalize();
+ subA2.normalize();
}
CanonicalForm
-reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
+reverseSubstReciproFq (const zz_pEX& F, const zz_pEX& G, int d, int k,
+ const Variable& alpha)
{
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);
-
+ zz_pEX f= F;
+ zz_pEX g= G;
+ int degf= deg(f);
+ int degg= deg(g);
- nmod_poly_t buf1,buf2, buf3;
+ zz_pEX buf1;
+ zz_pEX buf2;
+ zz_pEX buf3;
- if (nmod_poly_length (f) < (long) d*(k+1)) //zero padding
- nmod_poly_fit_length (f,(long)d*(k+1));
+ 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;
@@ -1398,6 +1378,8 @@ reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int 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)
@@ -1406,13 +1388,14 @@ reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
repLengthBuf1= 0;
else
repLengthBuf1= degfSubLf + 1;
- nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
+ buf1.rep.SetLength((long) repLengthBuf1);
+ buf1p= buf1.rep.elts();
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);
+ buf1p [ind]= fp [ind + lf];
+ buf1.normalize();
- repLengthBuf1= nmod_poly_length (buf1);
+ repLengthBuf1= buf1.rep.length();
if (deggSubLg >= d - 1)
repLengthBuf2= d - 1;
@@ -1421,23 +1404,27 @@ reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
else
repLengthBuf2= deggSubLg + 1;
- nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
+ buf2.rep.SetLength ((long) repLengthBuf2);
+ buf2p= buf2.rep.elts();
for (ind= 0; ind < repLengthBuf2; ind++)
- nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
+ buf2p [ind]= gp [ind + lg];
+ buf2.normalize();
- _nmod_poly_normalise (buf2);
- repLengthBuf2= nmod_poly_length (buf2);
+ repLengthBuf2= buf2.rep.length();
- nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
+ buf3.rep.SetLength((long) repLengthBuf2 + d);
+ buf3p= buf3.rep.elts();
+ buf2p= buf2.rep.elts();
+ buf1p= buf1.rep.elts();
for (ind= 0; ind < repLengthBuf1; ind++)
- nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
+ buf3p [ind]= buf1p [ind];
for (ind= repLengthBuf1; ind < d; ind++)
- nmod_poly_set_coeff_ui (buf3, ind, 0);
+ buf3p [ind]= zzpEZero;
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);
+ buf3p [ind + d]= buf2p [ind];
+ buf3.normalize();
- result += convertnmod_poly_t2FacCF (buf3, x)*power (y, i);
+ result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
i++;
@@ -1447,67 +1434,55 @@ reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
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++)
- 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()
- )
- );
+ gp [ind + lg] -= buf1p [ind];
}
+
if (lg < 0)
- {
- nmod_poly_clear (buf1);
- nmod_poly_clear (buf2);
- nmod_poly_clear (buf3);
break;
- }
+
+ buf2p= buf2.rep.elts();
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()
- )
- );
+ fp [ind + lf] -= buf2p [ind];
}
- 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)
+reverseSubstReciproFp (const zz_pX& F, const zz_pX& 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);
+ 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;
- fmpz_poly_t buf1,buf2, buf3;
+ zz_p *buf1p;
+ zz_p *buf2p;
+ zz_p *buf3p;
- if (fmpz_poly_length (f) < (long) d*(k+1)) //zero padding
- fmpz_poly_fit_length (f,(long)d*(k+1));
+ 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;
@@ -1515,7 +1490,7 @@ reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
int degfSubLf= degf;
int deggSubLg= degg-lg;
int repLengthBuf2, repLengthBuf1, ind, tmp;
- fmpz_t tmp1, tmp2;
+ zz_p zzpZero= zz_p();
while (degf >= lf || lg >= 0)
{
if (degfSubLf >= d)
@@ -1524,17 +1499,14 @@ reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
repLengthBuf1= 0;
else
repLengthBuf1= degfSubLf + 1;
+ buf1.rep.SetLength((long) repLengthBuf1);
- fmpz_poly_init2 (buf1, repLengthBuf1);
-
+ buf1p= buf1.rep.elts();
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);
+ buf1p [ind]= fp [ind + lf];
+ buf1.normalize();
- repLengthBuf1= fmpz_poly_length (buf1);
+ repLengthBuf1= buf1.rep.length();
if (deggSubLg >= d - 1)
repLengthBuf2= d - 1;
@@ -1543,33 +1515,29 @@ reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
else
repLengthBuf2= deggSubLg + 1;
- fmpz_poly_init2 (buf2, repLengthBuf2);
-
+ buf2.rep.SetLength ((long) repLengthBuf2);
+ buf2p= buf2.rep.elts();
for (ind= 0; ind < repLengthBuf2; ind++)
- {
- fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
- fmpz_poly_set_coeff_fmpz (buf2, ind, tmp1);
- }
+ buf2p [ind]= gp [ind + lg];
- _fmpz_poly_normalise (buf2);
- repLengthBuf2= fmpz_poly_length (buf2);
+ buf2.normalize();
- fmpz_poly_init2 (buf3, repLengthBuf2 + d);
+ 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++)
- {
- 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);
- }
+ buf3p [ind]= buf1p [ind];
for (ind= repLengthBuf1; ind < d; ind++)
- fmpz_poly_set_coeff_ui (buf3, ind, 0);
+ buf3p [ind]= zzpZero;
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);
+ buf3p [ind + d]= buf2p [ind];
+ buf3.normalize();
- result += convertFmpz_poly_t2FacCF (buf3, x)*power (y, i);
+ result += convertNTLzzpX2CF (buf3, x)*power (y, i);
i++;
@@ -1579,63 +1547,80 @@ reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
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++)
- {
- 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);
- }
+ gp [ind + lg] -= buf1p [ind];
}
if (lg < 0)
- {
- fmpz_poly_clear (buf1);
- fmpz_poly_clear (buf2);
- fmpz_poly_clear (buf3);
break;
- }
+
+ buf2p= buf2.rep.elts();
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);
- }
+ fp [ind + lf] -= buf2p [ind];
}
- 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 reverseSubstFq (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 nmod_poly_t F, int d)
+CanonicalForm reverseSubstFp (const zz_pX& 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);
+ zz_pX f= F;
+ zz_p *fp= f.rep.elts();
- nmod_poly_t buf;
+ zz_pX buf;
+ zz_p *bufp;
CanonicalForm result= 0;
int i= 0;
- int degf= nmod_poly_degree(f);
+ int degf= deg(f);
int k= 0;
int degfSubK, repLength, j;
while (degf >= k)
@@ -1646,67 +1631,56 @@ CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
else
repLength= degfSubK + 1;
- nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
+ buf.rep.SetLength ((long) repLength);
+ bufp= buf.rep.elts();
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);
+ bufp [j]= fp [j + k];
+ buf.normalize();
- result += convertnmod_poly_t2FacCF (buf, x)*power (y, i);
+ result += convertNTLzzpX2CF (buf, x)*power (y, i);
i++;
k= d*i;
- nmod_poly_clear (buf);
}
- nmod_poly_clear (f);
return result;
}
+// assumes input to be reduced mod M and to be an element of Fq not Fp
CanonicalForm
-mulMod2FLINTFpReci (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M)
+mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
+ CanonicalForm& M)
{
- int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
+ int d1= 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);
+ zz_pX F1, F2;
+ kronSubReciproFp (F1, F2, F, d1);
+ zz_pX G1, G2;
+ kronSubReciproFp (G1, G2, G, d1);
int k= d1*degree (M);
- nmod_poly_mullow (F1, F1, G1, (long) k);
+ MulTrunc (F1, F1, G1, (long) k);
- int degtailF= degree (tailcoeff (F), 1);;
+ 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);
- 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);
+ reverse (F2, F2);
+ reverse (G2, G2);
+ MulTrunc (F2, F2, G2, b + 1);
+ reverse (F2, F2, b);
- nmod_poly_clear (F1);
- nmod_poly_clear (F2);
- nmod_poly_clear (G1);
- nmod_poly_clear (G2);
- return result;
+ int d2= tmax (deg (F2)/d1, deg (F1)/d1);
+ return reverseSubstReciproFp (F1, F2, d1, d2);
}
+//Kronecker substitution
CanonicalForm
-mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M)
+mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
+ CanonicalForm& M)
{
CanonicalForm A= F;
CanonicalForm B= G;
@@ -1719,97 +1693,100 @@ mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
int d2= tmax (degAy, degBy);
if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
- return mulMod2FLINTFpReci (A, B, M);
+ return mulMod2NTLFpReci (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);
+ zz_pX NTLA= kronSubFp (A, d1);
+ zz_pX NTLB= kronSubFp (B, d1);
int k= d1*degree (M);
- nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
+ MulTrunc (NTLA, NTLA, NTLB, (long) k);
- A= reverseSubstFp (FLINTA, d1);
+ A= reverseSubstFp (NTLA, d1);
- nmod_poly_clear (FLINTA);
- nmod_poly_clear (FLINTB);
return A;
}
+// assumes input to be reduced mod M and to be an element of Fq not Fp
CanonicalForm
-mulMod2FLINTQReci (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M)
+mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
+ CanonicalForm& M, const Variable& alpha)
{
- int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
+ int d1= 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);
+ zz_pEX F1, F2;
+ kronSubReciproFq (F1, F2, F, d1, alpha);
+ zz_pEX G1, G2;
+ kronSubReciproFq (G1, G2, G, d1, alpha);
int k= d1*degree (M);
- fmpz_poly_mullow (F1, F1, G1, (long) k);
+ MulTrunc (F1, F1, G1, (long) k);
- int degtailF= degree (tailcoeff (F), 1);;
+ 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);
- 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);
+ reverse (F2, F2);
+ reverse (G2, G2);
+ MulTrunc (F2, F2, G2, b + 1);
+ reverse (F2, F2, b);
- fmpz_poly_clear (F1);
- fmpz_poly_clear (F2);
- fmpz_poly_clear (G1);
- fmpz_poly_clear (G2);
- return result;
+ int d2= tmax (deg (F2)/d1, deg (F1)/d1);
+ return reverseSubstReciproFq (F1, F2, d1, d2, alpha);
}
+#ifdef HAVE_FLINT
CanonicalForm
-mulMod2FLINTQ (const CanonicalForm& F, const CanonicalForm& G, const
- CanonicalForm& M)
+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;
- int degAx= degree (A, 1);
- int degBx= degree (B, 1);
- int d1= degAx + 1 + degBx;
+ 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);
- CanonicalForm f= bCommonDen (F);
- CanonicalForm g= bCommonDen (G);
- A *= f;
- B *= g;
+ int degMipo= degree (getMipo (alpha));
+ if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy) &&
+ (2*degAy > degree (M)))
+ return mulMod2NTLFqReci (A, B, M, alpha);
- 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);
+ zz_pEX NTLA= kronSubFq (A, d1, alpha);
+ zz_pEX NTLB= kronSubFq (B, d1, alpha);
- fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
- A= reverseSubstQ (FLINTA, d1);
- fmpz_poly_clear (FLINTA);
- fmpz_poly_clear (FLINTB);
- return A/(f*g);
-}
+ int k= d1*degree (M);
+
+ MulTrunc (NTLA, NTLA, NTLB, (long) k);
+ A= reverseSubstFq (NTLA, d1, alpha);
+
+ return A;
+ }
+ else
+#ifdef HAVE_FLINT
+ return mulMod2FLINTFp (A, B, M);
+#else
+ return mulMod2NTLFp (A, B, M);
#endif
+}
CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
const CanonicalForm& M)
@@ -1886,6 +1863,10 @@ CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
DEBOUTLN (cerr, "fatal end in mulMod2");
}
+// end bivariate polys
+//**********************
+// multivariate polys
+
CanonicalForm mod (const CanonicalForm& F, const CFList& M)
{
CanonicalForm A= F;
@@ -2030,6 +2011,10 @@ CanonicalForm prodMod (const CFList& L, const CFList& M)
}
}
+// end multivariate polys
+//***************************
+// division
+
CanonicalForm reverse (const CanonicalForm& F, int d)
{
if (d == 0)
@@ -2509,4 +2494,6 @@ void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
return;
}
+// end division
+
#endif
--
an open source computer algebra system
More information about the debian-science-commits
mailing list