[SCM] an open source computer algebra system branch, cleanedupstream, updated. 6125e540ca6d66c307958938a9d53b245507c323
Bernhard R. Link
brlink at debian.org
Tue Apr 24 15:54:44 UTC 2012
The following commit has been merged in the cleanedupstream branch:
commit 8f634ef6921d741300dffd1166bda24009c72f49
Author: Martin Lee <martinlee84 at web.de>
Date: Fri Jan 6 18:21:28 2012 +0100
chg: more changes in order to replace old factorization over integers
diff --git a/factory/facBivar.cc b/factory/facBivar.cc
index aeb25ec..f1d3d6b 100644
--- a/factory/facBivar.cc
+++ b/factory/facBivar.cc
@@ -524,6 +524,12 @@ CFList biFactorize (const CanonicalForm& F, const Variable& v)
return factors;
}
+ if (!extension)
+ {
+ for (CFListIterator i= uniFactors; i.hasItem(); i++)
+ i.getItem() /= lc (i.getItem());
+ }
+
A= A (y + evaluation, y);
int liftBound= degree (A, y) + 1;
@@ -531,14 +537,24 @@ CFList biFactorize (const CanonicalForm& F, const Variable& v)
bool earlySuccess= false;
CFList earlyFactors;
TIMING_START (fac_hensel_lift);
+ //out_cf ("A before= ",A, "\n");
uniFactors= henselLiftAndEarly0
(A, earlySuccess, earlyFactors, degs, liftBound,
uniFactors, evaluation);
TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
DEBOUTLN (cerr, "lifted factors= " << uniFactors);
+ //printf ("earlyFactors.length()= %d\n", earlyFactors.length());
+ //printf ("liftBound after= %d\n", liftBound);
+ //printf ("earlySuccess= %d\n", earlySuccess);
CanonicalForm MODl= power (y, liftBound);
+ CanonicalForm test= prod (uniFactors);
+ test *= LC (A,1);
+ test= mod (test, MODl);
+ //printf ("test == A %d\n", test == A);
+ //out_cf ("test= ", test, "\n");
+ //out_cf ("A= ", A, "\n");
factors= factorRecombination (uniFactors, A, MODl, degs, 1,
uniFactors.length()/2);
diff --git a/factory/facFactorize.cc b/factory/facFactorize.cc
index 1b24a7a..e75e731 100644
--- a/factory/facFactorize.cc
+++ b/factory/facFactorize.cc
@@ -599,7 +599,11 @@ multiFactorize (const CanonicalForm& F, const Variable& v)
//univariate case
if (F.isUnivariate())
{
- CFList result= conv (factorize (F, v));
+ CFList result;
+ if (v.level() != 1)
+ result= conv (factorize (F, v));
+ else
+ result= conv (factorize (F,true));
if (result.getFirst().inCoeffDomain())
result.removeFirst();
return result;
@@ -650,7 +654,10 @@ multiFactorize (const CanonicalForm& F, const Variable& v)
CFList factors;
if (A.isUnivariate ())
{
- factors= conv (factorize (A, v));
+ if (v.level() != 1)
+ factors= conv (factorize (A, v));
+ else
+ factors= conv (factorize (A,true));
append (factors, contentAFactors);
decompress (factors, N);
return factors;
diff --git a/factory/facHensel.cc b/factory/facHensel.cc
index 979a025..53ca5b7 100644
--- a/factory/facHensel.cc
+++ b/factory/facHensel.cc
@@ -559,10 +559,19 @@ mulNTL (const CanonicalForm& F, const CanonicalForm& G)
if ((!F.inCoeffDomain() && !G.inCoeffDomain()) && (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha)))
{
CanonicalForm result= mulFLINTQa (F, G, alpha);
+ CanonicalForm test= F*G;
+ if (test != result)
+ printf ("FAILLLLLLLL\n");
return result;
}
else if (!F.inCoeffDomain() && !G.inCoeffDomain())
+ {
+ CanonicalForm result= mulFLINTQ (F,G);
+ CanonicalForm test= F*G;
+ if (test != result)
+ printf ("FAILLLLLLLL2\n");
return mulFLINTQ (F, G);
+ }
#endif
return F*G;
}
@@ -615,7 +624,18 @@ modNTL (const CanonicalForm& F, const CanonicalForm& G)
if (getCharacteristic() == 0)
{
#ifdef HAVE_FLINT
- return modFLINTQ (F, G);
+ Variable alpha;
+ if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
+ {
+ CanonicalForm result= modFLINTQ (F,G);
+ CanonicalForm test= mod (F, G);
+ if (test != result)
+ printf ("FAILINMOD\n");
+ return result;
+ //return modFLINTQ (F, G);
+ }
+ else
+ return mod (F, G);
#else
return mod (F, G);
#endif
@@ -670,7 +690,18 @@ divNTL (const CanonicalForm& F, const CanonicalForm& G)
if (getCharacteristic() == 0)
{
#ifdef HAVE_FLINT
- return divFLINTQ (F,G);
+ Variable alpha;
+ if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
+ {
+ CanonicalForm result= divFLINTQ (F,G);
+ CanonicalForm test= div (F, G);
+ if (test != result)
+ printf ("FAILINDIV\n");
+ return result;
+ //return divFLINTQ (F,G);
+ }
+ else
+ return div (F, G);
#else
return div (F, G);
#endif
@@ -1407,6 +1438,12 @@ mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
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)
@@ -1444,9 +1481,647 @@ mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
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;
+
+ for (CFIterator i= A; i.hasTerms(); i++)
+ {
+ convertFacCF2nmod_poly_t (buf, i.coeff());
+
+ int k= i.exp()*d;
+ int kk= (degAy - i.exp())*d;
+ int bufRepLength= (int) nmod_poly_length (buf);
+ for (int 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())); //TODO muß man subA1 coeffs von subA1 vorher auf null setzen??
+ 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;
+
+ for (CFIterator i= A; i.hasTerms(); i++)
+ {
+ convertFacCF2Fmpz_poly_t (buf, i.coeff());
+
+ int k= i.exp()*d;
+ int kk= (degAy - i.exp())*d;
+ int bufRepLength= (int) fmpz_poly_length (buf);
+ for (int j= 0; j < bufRepLength; j++)
+ {
+ fmpz_poly_get_coeff_fmpz (coeff1, subA1, j+k); //TODO muß man subA1 coeffs von subA1 vorher auf null setzen??
+ // vielleicht ist es schneller fmpz_poly_get_ptr zu nehmen
+ 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;
+ int repLength;
+ 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 (int 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;
+ /*nmod_poly_init_preinv (buf1, getCharacteristic(), ninv);
+ nmod_poly_init_preinv (buf2, getCharacteristic(), ninv);
+ nmod_poly_init_preinv (buf3, getCharacteristic(), ninv);*/
+
+ 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;
+ int repLengthBuf1;
+ //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;
+ //cout << "repLengthBuf1= " << repLengthBuf1 << "\n";
+ //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
+ nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
+ //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
+
+ for (int ind= 0; ind < repLengthBuf1; ind++)
+ nmod_poly_set_coeff_ui (buf1, ind, nmod_poly_get_coeff_ui (f, ind+lf));
+ _nmod_poly_normalise (buf1);
+ //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
+
+ repLengthBuf1= nmod_poly_length (buf1);
+
+ if (deggSubLg >= d - 1)
+ repLengthBuf2= d - 1;
+ else if (deggSubLg < 0)
+ repLengthBuf2= 0;
+ else
+ repLengthBuf2= deggSubLg + 1;
+
+ //cout << "repLengthBuf2= " << repLengthBuf2 << "\n";
+ //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
+ nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
+ //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
+ for (int ind= 0; ind < repLengthBuf2; ind++)
+ nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
+
+ _nmod_poly_normalise (buf2);
+ //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
+ repLengthBuf2= nmod_poly_length (buf2);
+
+ nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
+ for (int ind= 0; ind < repLengthBuf1; ind++)
+ nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
+ for (int ind= repLengthBuf1; ind < d; ind++)
+ nmod_poly_set_coeff_ui (buf3, ind, 0);
+ for (int 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;
+ int tmp= tmin (repLengthBuf1, deggSubLg + 1);
+ for (int 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()));
+ //cout << "nmod_poly_get_coeff_ui (g, ind + lg)= " << nmod_poly_get_coeff_ui (g, ind + lg) << "\n";
+ //}
+ }
+ if (lg < 0)
+ {
+ nmod_poly_clear (buf1);
+ nmod_poly_clear (buf2);
+ nmod_poly_clear (buf3);
+ break;
+ }
+ if (degfSubLf >= 0)
+ {
+ for (int 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_init_preinv (buf1, getCharacteristic(), ninv);
+ nmod_poly_init_preinv (buf2, getCharacteristic(), ninv);
+ nmod_poly_init_preinv (buf3, getCharacteristic(), ninv);*/
+ }
+
+ /*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;
+ /*nmod_poly_init_preinv (buf1, getCharacteristic(), ninv);
+ nmod_poly_init_preinv (buf2, getCharacteristic(), ninv);
+ nmod_poly_init_preinv (buf3, getCharacteristic(), ninv);*/
+
+ 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;
+ int repLengthBuf1;
+ //zz_p zzpZero= zz_p();
+ fmpz_t tmp1, tmp2;
+ while (degf >= lf || lg >= 0)
+ {
+ if (degfSubLf >= d)
+ repLengthBuf1= d;
+ else if (degfSubLf < 0)
+ repLengthBuf1= 0;
+ else
+ repLengthBuf1= degfSubLf + 1;
+ //cout << "repLengthBuf1= " << repLengthBuf1 << "\n";
+ //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
+ fmpz_poly_init2 (buf1, repLengthBuf1);
+ //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
+
+ for (int 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);
+ //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
+
+ repLengthBuf1= fmpz_poly_length (buf1);
+
+ if (deggSubLg >= d - 1)
+ repLengthBuf2= d - 1;
+ else if (deggSubLg < 0)
+ repLengthBuf2= 0;
+ else
+ repLengthBuf2= deggSubLg + 1;
+
+ //cout << "repLengthBuf2= " << repLengthBuf2 << "\n";
+ //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
+ fmpz_poly_init2 (buf2, repLengthBuf2);
+ //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
+ for (int 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);
+ //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
+ repLengthBuf2= fmpz_poly_length (buf2);
+
+ fmpz_poly_init2 (buf3, repLengthBuf2 + d);
+ for (int 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 (int ind= repLengthBuf1; ind < d; ind++)
+ fmpz_poly_set_coeff_ui (buf3, ind, 0);
+ for (int 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;
+ int tmp= tmin (repLengthBuf1, deggSubLg + 1);
+ for (int 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);
+ //cout << "nmod_poly_get_coeff_ui (g, ind + lg)= " << nmod_poly_get_coeff_ui (g, ind + lg) << "\n";
+ }
+ }
+ if (lg < 0)
+ {
+ fmpz_poly_clear (buf1);
+ fmpz_poly_clear (buf2);
+ fmpz_poly_clear (buf3);
+ break;
+ }
+ if (degfSubLf >= 0)
+ {
+ for (int 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);
+ /*nmod_poly_init_preinv (buf1, getCharacteristic(), ninv);
+ nmod_poly_init_preinv (buf2, getCharacteristic(), ninv);
+ nmod_poly_init_preinv (buf3, getCharacteristic(), ninv);*/
+ }
+
+ /*nmod_poly_clear (buf1);
+ nmod_poly_clear (buf2);
+ nmod_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;
+ //nmod_poly_init_preinv (buf, getCharacteristic(), ninv);
+ CanonicalForm result= 0;
+ int i= 0;
+ int degf= nmod_poly_degree(f);
+ int k= 0;
+ int degfSubK;
+ int repLength;
+ while (degf >= k)
+ {
+ degfSubK= degf - k;
+ if (degfSubK >= d)
+ repLength= d;
+ else
+ repLength= degfSubK + 1;
+
+ nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
+ for (int j= 0; j < repLength; j++)
+ nmod_poly_set_coeff_ui (buf, j, nmod_poly_get_coeff_ui (f, j + k));
+ //printf ("after set\n");
+ _nmod_poly_normalise (buf);
+ //printf ("after normalise\n");
+
+ result += convertnmod_poly_t2FacCF (buf, x)*power (y, i);
+ i++;
+ k= d*i;
+ nmod_poly_clear (buf);
+ //nmod_poly_init_preinv (buf, getCharacteristic(), ninv);
+ }
+ //printf ("after while\n");
+ //nmod_poly_clear (buf);
+ nmod_poly_clear (f);
+ //printf ("after clear\n");
+
+ 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);
+ /*zz_pX TF1, TF2;
+ kronSubRecipro (TF1, TF2, F, d1);
+ Variable x= Variable (1);
+ cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF1, x) == convertnmod_poly_t2FacCF (F1, x)) << "\n";
+ cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF2, x) == convertnmod_poly_t2FacCF (F2, x)) << "\n";*/
+
+ 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 TG1, TG2;
+ kronSubRecipro (TG1, TG2, G, d1);
+ cout << "testkronSubReci= " << (convertNTLzzpX2CF (TG1, x) == convertnmod_poly_t2FacCF (G1, x)) << "\n";
+ cout << "testkronSubReci= " << (convertNTLzzpX2CF (TG2, x) == convertnmod_poly_t2FacCF (G2, x)) << "\n";*/
+
+
+ int k= d1*degree (M);
+ nmod_poly_mullow (F1, F1, G1, (long) k);
+ /*MulTrunc (TF1, TF1, TG1, (long) k);
+ cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF1, x) == convertnmod_poly_t2FacCF (F1, x)) << "\n";*/
+
+ 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);
+ /*mul (TF2, TF2, TG2);
+ if (deg (TF2) > k)
+ {
+ int b= deg (TF2) - k + 2;
+ TF2 >>= b;
+ }
+ cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF2, x) == convertnmod_poly_t2FacCF (F2, x)) << "\n";*/
+ int d2= tmax (nmod_poly_degree (F2)/d1, nmod_poly_degree (F1)/d1);
+
+
+ CanonicalForm result= reverseSubst (F1, F2, d1, d2);
+ //CanonicalForm NTLres= reverseSubst (TF1, TF2, d1, d2);
+ //cout << "FINALtestkronSubReci= " << (NTLres == result) << "\n";
+
+ 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);
+ /*zz_pX TF1, TF2;
+ kronSubRecipro (TF1, TF2, F, d1);
+ Variable x= Variable (1);
+ cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF1, x) == convertnmod_poly_t2FacCF (F1, x)) << "\n";
+ cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF2, x) == convertnmod_poly_t2FacCF (F2, x)) << "\n";*/
+
+ fmpz_poly_t G1, G2;
+ fmpz_poly_init (G1);
+ fmpz_poly_init (G2);
+ kronSubRecipro (G1, G2, G, d1);
+ /*zz_pX TG1, TG2;
+ kronSubRecipro (TG1, TG2, G, d1);
+ cout << "testkronSubReci= " << (convertNTLzzpX2CF (TG1, x) == convertnmod_poly_t2FacCF (G1, x)) << "\n";
+ cout << "testkronSubReci= " << (convertNTLzzpX2CF (TG2, x) == convertnmod_poly_t2FacCF (G2, x)) << "\n";*/
+
+
+ int k= d1*degree (M);
+ fmpz_poly_mullow (F1, F1, G1, (long) k);
+ /*MulTrunc (TF1, TF1, TG1, (long) k);
+ cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF1, x) == convertnmod_poly_t2FacCF (F1, x)) << "\n";*/
+
+ 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);
+ /*mul (TF2, TF2, TG2);
+ if (deg (TF2) > k)
+ {
+ int b= deg (TF2) - k + 2;
+ TF2 >>= b;
+ }
+ cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF2, x) == convertnmod_poly_t2FacCF (F2, x)) << "\n";*/
+ int d2= tmax (fmpz_poly_degree (F2)/d1, fmpz_poly_degree (F1)/d1);
+
+
+ CanonicalForm result= reverseSubst (F1, F2, d1, d2);
+
+ //CanonicalForm NTLres= reverseSubst (TF1, TF2, d1, d2);
+ //cout << "FINALtestkronSubReci= " << (NTLres == result) << "\n";
+
+ 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 degAy= degree (A, 2);
+ int degBx= degree (B, 1);
+ //int degBy= degree (B, 2);
+ int d1= degAx + 1 + degBx;
+ //int d2= tmax (degAy, degBy);
+
+ 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)
{
@@ -1482,6 +2157,19 @@ CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
if (sizeF < fallBackToNaive || sizeG < fallBackToNaive)
return mod (F*G, M);
+#ifdef HAVE_FLINT
+ Variable alpha;
+ if (getCharacteristic() == 0 && !hasFirstAlgVar (F, alpha) && ! hasFirstAlgVar (G, alpha))
+ {
+ CanonicalForm result= mulMod2FLINTQ (F, G, M);
+ CanonicalForm test= mod (F*G, M);
+ if (test != result)
+ printf ("GOTTVERDAMMT\n");
+ return result;
+ //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);
--
an open source computer algebra system
More information about the debian-science-commits
mailing list