[SCM] an open source computer algebra system branch, cleanedupstream, updated. 6125e540ca6d66c307958938a9d53b245507c323
Bernhard R. Link
brlink at debian.org
Tue Apr 24 15:54:42 UTC 2012
The following commit has been merged in the cleanedupstream branch:
commit a5820654ccaf1bc2ef713447d94c147765e9f8ac
Author: Martin Lee <martinlee84 at web.de>
Date: Tue Dec 20 17:08:05 2011 +0000
chg: switched off old factory factorization over Z
chg: added some function declarations to FLINTconvert.h
chg: added a lot of modular multiplication code and univariate
arithmetic over Q and Q(a) using FLINT
diff --git a/factory/FLINTconvert.h b/factory/FLINTconvert.h
index 8058fbf..09d5fee 100644
--- a/factory/FLINTconvert.h
+++ b/factory/FLINTconvert.h
@@ -12,6 +12,7 @@ extern "C"
#include <fmpz.h>
#include <fmpq.h>
#include <fmpz_poly.h>
+#include <fmpq_poly.h>
#include <nmod_poly.h>
#ifdef __cplusplus
}
@@ -24,6 +25,8 @@ CanonicalForm convertFmpz_poly_t2FacCF (fmpz_poly_t poly, const Variable& x);
void convertFacCF2nmod_poly_t (nmod_poly_t result, const CanonicalForm& f);
CanonicalForm convertnmod_poly_t2FacCF (nmod_poly_t poly, const Variable& x);
void convertCF2Fmpq (fmpq_t result, const CanonicalForm& f);
+CanonicalForm convertFmpq_poly_t2FacCF (fmpq_poly_t p, const Variable& x);
+void convertFacCF2Fmpq_poly_t (fmpq_poly_t result, const CanonicalForm& f);
CFFList convertFLINTnmod_poly_factor2FacCFFList (nmod_poly_factor_t fac,
mp_limb_t leadingCoeff,
const Variable& x
diff --git a/factory/cf_factor.cc b/factory/cf_factor.cc
index e2b813a..b1350aa 100644
--- a/factory/cf_factor.cc
+++ b/factory/cf_factor.cc
@@ -637,7 +637,17 @@ CFFList factorize ( const CanonicalForm & f, bool issqrfree )
}
else
{
- F = ZFactorizeMultivariate( fz, issqrfree );
+ On (SW_RATIONAL);
+ if (issqrfree)
+ {
+ CFList factors;
+ factors= ratSqrfFactorize (fz, Variable (1));
+ for (CFListIterator i= factors; i.hasItem(); i++)
+ F.append (CFFactor (i.getItem(), 1));
+ }
+ else
+ F = ratFactorize (fz, Variable (1));
+ Off (SW_RATIONAL);
}
if ( on_rational )
diff --git a/factory/facFactorize.cc b/factory/facFactorize.cc
index 632062b..1b24a7a 100644
--- a/factory/facFactorize.cc
+++ b/factory/facFactorize.cc
@@ -873,7 +873,7 @@ multiFactorize (const CanonicalForm& F, const Variable& v)
A /= hh;
- if (LucksWangSparseHeuristic (A, biFactors, 2, leadingCoeffs2 [A.level() - 3],
+ /*if (LucksWangSparseHeuristic (A, biFactors, 2, leadingCoeffs2 [A.level() - 3],
factors))
{
int check= factors.length();
@@ -888,7 +888,7 @@ multiFactorize (const CanonicalForm& F, const Variable& v)
}
else
factors= CFList();
- }
+ }*/
//shifting to zero
diff --git a/factory/facHensel.cc b/factory/facHensel.cc
index 55b13a7..979a025 100644
--- a/factory/facHensel.cc
+++ b/factory/facHensel.cc
@@ -31,6 +31,228 @@
#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);
+
+
+ /*fmpz_t coeff;
+ for (CFIterator i= A; i.hasTerms(); i++)
+ {
+ if (i.coeff().inBaseDomain())
+ {
+ convertCF2Fmpz (coeff, i.coeff());
+ fmpz_poly_set_coeff_fmpz (result, i.exp()*d, coeff);
+ }
+ else
+ {
+ for (CFIterator j= i.coeff();j.hasTerms(); j++)
+ {
+ convertCF2Fmpz (coeff, j.coeff());
+ fmpz_poly_set_coeff_fmpz (result, i.exp()*d+j.exp(), coeff);
+ }
+ }
+ }
+ fmpz_clear (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;
+ CanonicalForm coeff;
+ fmpz* tmp;
+ while (degf >= k)
+ {
+ coeff= 0;
+ degfSubK= degf - k;
+ if (degfSubK >= d)
+ repLength= d;
+ else
+ repLength= degfSubK + 1;
+
+ for (int 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;
+}
+
+#endif
+
static
CFList productsNTL (const CFList& factors, const CanonicalForm& M)
{
@@ -331,7 +553,19 @@ 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)
@@ -350,10 +584,20 @@ mulNTL (const CanonicalForm& F, const CanonicalForm& G)
}
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;
}
@@ -369,7 +613,13 @@ modNTL (const CanonicalForm& F, const CanonicalForm& G)
return mod (F,G);
if (getCharacteristic() == 0)
+ {
+#ifdef HAVE_FLINT
+ return modFLINTQ (F, G);
+#else
return mod (F, G);
+#endif
+ }
ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
ASSERT (F.level() == G.level(), "expected polys of same level");
@@ -389,10 +639,20 @@ modNTL (const CanonicalForm& F, const CanonicalForm& G)
}
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;
}
@@ -408,7 +668,13 @@ divNTL (const CanonicalForm& F, const CanonicalForm& G)
return div (F,G);
if (getCharacteristic() == 0)
+ {
+#ifdef HAVE_FLINT
+ return divFLINTQ (F,G);
+#else
return div (F, G);
+#endif
+ }
ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
ASSERT (F.level() == G.level(), "expected polys of same level");
@@ -428,10 +694,20 @@ divNTL (const CanonicalForm& F, const CanonicalForm& G)
}
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;
}
@@ -501,6 +777,101 @@ CanonicalForm mod (const CanonicalForm& F, const CFList& M)
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;
+
+ for (CFIterator i= A; i.hasTerms(); i++)
+ {
+ convertFacCF2nmod_poly_t (buf, i.coeff());
+
+ int k= i.exp()*d;
+ int bufRepLength= (int) nmod_poly_length (buf);
+ for (int 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);
--
an open source computer algebra system
More information about the debian-science-commits
mailing list