Skip to content

Commit 1d0c92d

Browse files
mmkleemmklee
committed
chg: use convex dense compression in gcd
1 parent 8243b8d commit 1d0c92d

File tree

3 files changed

+212
-9
lines changed

3 files changed

+212
-9
lines changed

‎factory/cfNewtonPolygon.cc‎

Lines changed: 14 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -594,11 +594,15 @@ void convexDense(int** points, int sizePoints, mat_ZZ& M, vec_ZZ& A)
594594
}
595595

596596
CanonicalForm
597-
compress (const CanonicalForm& F, mat_ZZ& M, vec_ZZ& A)
597+
compress (const CanonicalForm& F, mat_ZZ& M, vec_ZZ& A, bool computeMA)
598598
{
599599
int n;
600-
int ** newtonPolyg= newtonPolygon (F, n);
601-
convexDense (newtonPolyg, n, M, A);
600+
int ** newtonPolyg;
601+
if (computeMA)
602+
{
603+
newtonPolyg= newtonPolygon (F, n);
604+
convexDense (newtonPolyg, n, M, A);
605+
}
602606
CanonicalForm result= 0;
603607
ZZ expX, expY;
604608
Variable x= Variable (1);
@@ -675,11 +679,14 @@ compress (const CanonicalForm& F, mat_ZZ& M, vec_ZZ& A)
675679
result += Lc (tmp)*power (x, d);
676680
}
677681

678-
for (int i= 0; i < n; i++)
679-
delete [] newtonPolyg [i];
680-
delete [] newtonPolyg;
682+
if (computeMA)
683+
{
684+
for (int i= 0; i < n; i++)
685+
delete [] newtonPolyg [i];
686+
delete [] newtonPolyg;
687+
M= inv (M);
688+
}
681689

682-
M= inv (M);
683690
return result;
684691
}
685692

‎factory/cfNewtonPolygon.h‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -81,8 +81,10 @@ void convexDense (int** points, ///< [in, out] a set of points in Z^2, returns
8181
CanonicalForm
8282
compress (const CanonicalForm& F, ///< [in] compressed, i.e. F.level()==2,
8383
///< bivariate poly
84-
mat_ZZ& inverseM, ///< [in,out] returns the inverse of M
85-
vec_ZZ& A ///< [in,out] returns translation
84+
mat_ZZ& inverseM, ///< [in,out] returns the inverse of M,
85+
///< if computeMA==true, M otherwise
86+
vec_ZZ& A, ///< [in,out] returns translation
87+
bool computeMA= true ///< [in] whether to compute M and A
8688
);
8789

8890
/// decompress a bivariate poly

‎factory/cf_gcd_smallp.cc‎

Lines changed: 194 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@
2929
#include "cf_random.h"
3030
#include "cf_reval.h"
3131
#include "facHensel.h"
32+
#include "cfNewtonPolygon.h"
3233

3334
// iinline helper functions:
3435
#include "cf_map_ext.h"
@@ -593,6 +594,67 @@ GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G,
593594
ppB= B/cB;
594595
}
595596

597+
int sizeNewtonPolyg;
598+
int ** newtonPolyg= NULL;
599+
mat_ZZ MM;
600+
vec_ZZ V;
601+
bool compressConvexDense= (ppA.level() == 2 && ppB.level() == 2);
602+
if (compressConvexDense)
603+
{
604+
cA= content (ppA, 1);
605+
cB= content (ppB, 1);
606+
ppA /= cA;
607+
ppB /= cB;
608+
gcdcAcB *= gcd (cA, cB);
609+
if (ppA.isUnivariate() || ppB.isUnivariate())
610+
{
611+
if (ppA.level() == ppB.level())
612+
{
613+
if (substitute > 1)
614+
return N (reverseSubst (gcd (ppA, ppB)*gcdcAcB, substitute));
615+
else
616+
return N (gcd (ppA, ppB)*gcdcAcB);
617+
}
618+
else
619+
{
620+
if (substitute > 1)
621+
return N (reverseSubst (gcdcAcB, substitute));
622+
else
623+
return N (gcdcAcB);
624+
}
625+
}
626+
627+
newtonPolyg= newtonPolygon (ppA,ppB, sizeNewtonPolyg);
628+
convexDense (newtonPolyg, sizeNewtonPolyg, MM, V);
629+
630+
for (int i= 0; i < sizeNewtonPolyg; i++)
631+
delete [] newtonPolyg[i];
632+
delete [] newtonPolyg;
633+
634+
ppA= compress (ppA, MM, V, false);
635+
ppB= compress (ppB, MM, V, false);
636+
MM= inv (MM);
637+
638+
if (ppA.isUnivariate() && ppB.isUnivariate())
639+
{
640+
if (ppA.level() == ppB.level())
641+
{
642+
if (substitute > 1)
643+
return N (reverseSubst (decompress (gcd (ppA, ppB), MM, V)*gcdcAcB,
644+
substitute));
645+
else
646+
return N (decompress (gcd (ppA, ppB), MM, V)*gcdcAcB);
647+
}
648+
else
649+
{
650+
if (substitute > 1)
651+
return N (reverseSubst (gcdcAcB, substitute));
652+
else
653+
return N (gcdcAcB);
654+
}
655+
}
656+
}
657+
596658
CanonicalForm lcA, lcB; // leading coefficients of A and B
597659
CanonicalForm gcdlcAlcB;
598660

@@ -778,6 +840,8 @@ GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G,
778840
ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
779841
ppH /= Lc(ppH);
780842
DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
843+
if (compressConvexDense)
844+
ppH= decompress (ppH, MM, V);
781845
if (substitute > 1)
782846
{
783847
ppH= reverseSubst (ppH, substitute);
@@ -788,6 +852,8 @@ GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G,
788852
}
789853
else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
790854
{
855+
if (compressConvexDense)
856+
ppH= decompress (ppH, MM, V);
791857
if (substitute > 1)
792858
{
793859
ppH= reverseSubst (ppH, substitute);
@@ -900,6 +966,67 @@ CanonicalForm GCD_GF (const CanonicalForm& F, const CanonicalForm& G,
900966
ppB= B/cB;
901967
}
902968

969+
int sizeNewtonPolyg;
970+
int ** newtonPolyg= NULL;
971+
mat_ZZ MM;
972+
vec_ZZ V;
973+
bool compressConvexDense= (ppA.level() == 2 && ppB.level() == 2);
974+
if (compressConvexDense)
975+
{
976+
cA= content (ppA, 1);
977+
cB= content (ppB, 1);
978+
ppA /= cA;
979+
ppB /= cB;
980+
gcdcAcB *= gcd (cA, cB);
981+
if (ppA.isUnivariate() || ppB.isUnivariate())
982+
{
983+
if (ppA.level() == ppB.level())
984+
{
985+
if (substitute > 1)
986+
return N (reverseSubst (gcd (ppA, ppB)*gcdcAcB, substitute));
987+
else
988+
return N (gcd (ppA, ppB)*gcdcAcB);
989+
}
990+
else
991+
{
992+
if (substitute > 1)
993+
return N (reverseSubst (gcdcAcB, substitute));
994+
else
995+
return N (gcdcAcB);
996+
}
997+
}
998+
999+
newtonPolyg= newtonPolygon (ppA,ppB, sizeNewtonPolyg);
1000+
convexDense (newtonPolyg, sizeNewtonPolyg, MM, V);
1001+
1002+
for (int i= 0; i < sizeNewtonPolyg; i++)
1003+
delete [] newtonPolyg[i];
1004+
delete [] newtonPolyg;
1005+
1006+
ppA= compress (ppA, MM, V, false);
1007+
ppB= compress (ppB, MM, V, false);
1008+
MM= inv (MM);
1009+
1010+
if (ppA.isUnivariate() && ppB.isUnivariate())
1011+
{
1012+
if (ppA.level() == ppB.level())
1013+
{
1014+
if (substitute > 1)
1015+
return N (reverseSubst (decompress (gcd (ppA, ppB), MM, V)*gcdcAcB,
1016+
substitute));
1017+
else
1018+
return N (decompress (gcd (ppA, ppB), MM, V)*gcdcAcB);
1019+
}
1020+
else
1021+
{
1022+
if (substitute > 1)
1023+
return N (reverseSubst (gcdcAcB, substitute));
1024+
else
1025+
return N (gcdcAcB);
1026+
}
1027+
}
1028+
}
1029+
9031030
CanonicalForm lcA, lcB; // leading coefficients of A and B
9041031
CanonicalForm gcdlcAlcB;
9051032

@@ -1066,6 +1193,8 @@ CanonicalForm GCD_GF (const CanonicalForm& F, const CanonicalForm& G,
10661193
DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
10671194
ppH= GFMapDown (ppH, k);
10681195
DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
1196+
if (compressConvexDense)
1197+
ppH= decompress (ppH, MM, V);
10691198
if (substitute > 1)
10701199
{
10711200
ppH= reverseSubst (ppH, substitute);
@@ -1079,6 +1208,8 @@ CanonicalForm GCD_GF (const CanonicalForm& F, const CanonicalForm& G,
10791208
{
10801209
if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
10811210
{
1211+
if (compressConvexDense)
1212+
ppH= decompress (ppH, MM, V);
10821213
if (substitute > 1)
10831214
{
10841215
ppH= reverseSubst (ppH, substitute);
@@ -1199,6 +1330,67 @@ CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm& G,
11991330
ppB= B/cB;
12001331
}
12011332

1333+
int sizeNewtonPolyg;
1334+
int ** newtonPolyg= NULL;
1335+
mat_ZZ MM;
1336+
vec_ZZ V;
1337+
bool compressConvexDense= (ppA.level() == 2 && ppB.level() == 2);
1338+
if (compressConvexDense)
1339+
{
1340+
cA= content (ppA, 1);
1341+
cB= content (ppB, 1);
1342+
ppA /= cA;
1343+
ppB /= cB;
1344+
gcdcAcB *= gcd (cA, cB);
1345+
if (ppA.isUnivariate() || ppB.isUnivariate())
1346+
{
1347+
if (ppA.level() == ppB.level())
1348+
{
1349+
if (substitute > 1)
1350+
return N (reverseSubst (gcd (ppA, ppB)*gcdcAcB, substitute));
1351+
else
1352+
return N (gcd (ppA, ppB)*gcdcAcB);
1353+
}
1354+
else
1355+
{
1356+
if (substitute > 1)
1357+
return N (reverseSubst (gcdcAcB, substitute));
1358+
else
1359+
return N (gcdcAcB);
1360+
}
1361+
}
1362+
1363+
newtonPolyg= newtonPolygon (ppA,ppB, sizeNewtonPolyg);
1364+
convexDense (newtonPolyg, sizeNewtonPolyg, MM, V);
1365+
1366+
for (int i= 0; i < sizeNewtonPolyg; i++)
1367+
delete [] newtonPolyg[i];
1368+
delete [] newtonPolyg;
1369+
1370+
ppA= compress (ppA, MM, V, false);
1371+
ppB= compress (ppB, MM, V, false);
1372+
MM= inv (MM);
1373+
1374+
if (ppA.isUnivariate() && ppB.isUnivariate())
1375+
{
1376+
if (ppA.level() == ppB.level())
1377+
{
1378+
if (substitute > 1)
1379+
return N (reverseSubst (decompress (gcd (ppA, ppB), MM, V)*gcdcAcB,
1380+
substitute));
1381+
else
1382+
return N (decompress (gcd (ppA, ppB), MM, V)*gcdcAcB);
1383+
}
1384+
else
1385+
{
1386+
if (substitute > 1)
1387+
return N (reverseSubst (gcdcAcB, substitute));
1388+
else
1389+
return N (gcdcAcB);
1390+
}
1391+
}
1392+
}
1393+
12021394
CanonicalForm lcA, lcB; // leading coefficients of A and B
12031395
CanonicalForm gcdlcAlcB;
12041396
lcA= uni_lcoeff (ppA);
@@ -1418,6 +1610,8 @@ CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm& G,
14181610
DEBOUTLN (cerr, "ppH= " << ppH);
14191611
if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
14201612
{
1613+
if (compressConvexDense)
1614+
ppH= decompress (ppH, MM, V);
14211615
if (substitute > 1)
14221616
{
14231617
ppH= reverseSubst (ppH, substitute);

0 commit comments

Comments
 (0)