Skip to content

Commit

Permalink
Fix modular GCD for finite fields of small cardinality:
Browse files Browse the repository at this point in the history 
 - overcome the issue with LinZip #7
 - make switch between different x_n in `F[x_n][x_1, ..., x_(n-1)]` representation
  • Loading branch information
@PoslavskySV
PoslavskySV committed Sep 5, 2017
1 parent ff6d7ba commit 9a0d819
Show file tree
Hide file tree
Showing 2 changed files with 132 additions and 64 deletions.
167 changes: 104 additions & 63 deletions src/main/java/cc/r2/core/poly/multivar/MultivariateGCD.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -284,6 +284,8 @@ Poly trivialGCD(Poly a, Poly b) {
return gcdWithMonomial(a.lt(), b);
if (b.size() == 1)
return gcdWithMonomial(b.lt(), a);
if (a.equals(b))
return a.clone();
return null;
}

Expand Down Expand Up @@ -951,44 +953,56 @@ private static <E> MultivariatePolynomial<E> ModularGCDFiniteField(GCDInput<Mono
MultivariatePolynomial<E> a = gcdInput.aReduced;
MultivariatePolynomial<E> b = gcdInput.bReduced;

int uVariable = a.nVariables - 1;
if (a.domain instanceof IntegersModulo && a.coefficientDomainCardinality().isLong()) {
// use machine integers
@SuppressWarnings("unchecked")
MultivariatePolynomial<lUnivariatePolynomialZp>
ua = asOverUnivariate0((MultivariatePolynomial<BigInteger>) a, uVariable),
ub = asOverUnivariate0((MultivariatePolynomial<BigInteger>) b, uVariable);

lUnivariatePolynomialZp aContent = ua.content(), bContent = ub.content();
lUnivariatePolynomialZp contentGCD = ua.domain.gcd(aContent, bContent);

ua = ua.divideOrNull(aContent);
ub = ub.divideOrNull(bContent);

MultivariatePolynomial<lUnivariatePolynomialZp> ugcd =
ModularGCDFiniteField0(ua, ub, gcdInput.degreeBounds[uVariable], gcdInput.finiteExtensionDegree);
ugcd = ugcd.multiply(contentGCD);
// MultivariatePolynomial<E> pContentGCD = contentGCD(a, b, 0, MultivariateGCD::ModularGCDFiniteField);
// if (!pContentGCD.isConstant()) {
// a = divideExact(a, pContentGCD);
// b = divideExact(b, pContentGCD);
// return gcdInput.restoreGCD(ModularGCDFiniteField(a, b).multiply(pContentGCD));
// }

for (int uVariable = a.nVariables - 1; uVariable >= 0; --uVariable)
if (a.domain instanceof IntegersModulo && a.coefficientDomainCardinality().isLong()) {
// use machine integers
@SuppressWarnings("unchecked")
MultivariatePolynomial<lUnivariatePolynomialZp>
ua = asOverUnivariate0((MultivariatePolynomial<BigInteger>) a, uVariable),
ub = asOverUnivariate0((MultivariatePolynomial<BigInteger>) b, uVariable);

lUnivariatePolynomialZp aContent = ua.content(), bContent = ub.content();
lUnivariatePolynomialZp contentGCD = ua.domain.gcd(aContent, bContent);

ua = ua.divideOrNull(aContent);
ub = ub.divideOrNull(bContent);

MultivariatePolynomial<lUnivariatePolynomialZp> ugcd =
ModularGCDFiniteField0(ua, ub, gcdInput.degreeBounds[uVariable], gcdInput.finiteExtensionDegree);
if (ugcd == null)
continue;
ugcd = ugcd.multiply(contentGCD);

@SuppressWarnings("unchecked")
MultivariatePolynomial<E> r = gcdInput.restoreGCD((MultivariatePolynomial<E>) asNormalMultivariate0(ugcd, uVariable));
return r;
} else {
MultivariatePolynomial<UnivariatePolynomial<E>>
ua = a.asOverUnivariateEliminate(uVariable),
ub = b.asOverUnivariateEliminate(uVariable);
@SuppressWarnings("unchecked")
MultivariatePolynomial<E> r = gcdInput.restoreGCD((MultivariatePolynomial<E>) asNormalMultivariate0(ugcd, uVariable));
return r;
} else {
MultivariatePolynomial<UnivariatePolynomial<E>>
ua = a.asOverUnivariateEliminate(uVariable),
ub = b.asOverUnivariateEliminate(uVariable);

UnivariatePolynomial<E> aContent = ua.content(), bContent = ub.content();
UnivariatePolynomial<E> contentGCD = ua.domain.gcd(aContent, bContent);
UnivariatePolynomial<E> aContent = ua.content(), bContent = ub.content();
UnivariatePolynomial<E> contentGCD = ua.domain.gcd(aContent, bContent);

ua = ua.divideOrNull(aContent);
ub = ub.divideOrNull(bContent);
ua = ua.divideOrNull(aContent);
ub = ub.divideOrNull(bContent);

MultivariatePolynomial<UnivariatePolynomial<E>> ugcd =
ModularGCDFiniteField0(ua, ub, gcdInput.degreeBounds[uVariable], gcdInput.finiteExtensionDegree);
ugcd = ugcd.multiply(contentGCD);
MultivariatePolynomial<UnivariatePolynomial<E>> ugcd =
ModularGCDFiniteField0(ua, ub, gcdInput.degreeBounds[uVariable], gcdInput.finiteExtensionDegree);
if (ugcd == null)
continue;
ugcd = ugcd.multiply(contentGCD);

return gcdInput.restoreGCD(MultivariatePolynomial.asNormalMultivariate(ugcd, uVariable));
}
return gcdInput.restoreGCD(MultivariatePolynomial.asNormalMultivariate(ugcd, uVariable));
}
throw new RuntimeException();
}

private static MultivariatePolynomial<lUnivariatePolynomialZp> asOverUnivariate0(MultivariatePolynomial<BigInteger> poly, int variable) {
Expand Down Expand Up @@ -1054,25 +1068,37 @@ private static lMultivariatePolynomialZp lModularGCDFiniteField(GCDInput<lMonomi
lMultivariatePolynomialZp a = gcdInput.aReduced;
lMultivariatePolynomialZp b = gcdInput.bReduced;

int uVariable = a.nVariables - 1;
MultivariatePolynomial<lUnivariatePolynomialZp>
ua = a.asOverUnivariateEliminate(uVariable),
ub = b.asOverUnivariateEliminate(uVariable);
// lMultivariatePolynomialZp pContentGCD = contentGCD(a, b, 0, MultivariateGCD::ModularGCDFiniteField);
// if (!pContentGCD.isConstant()) {
// a = divideExact(a, pContentGCD);
// b = divideExact(b, pContentGCD);
// return gcdInput.restoreGCD(ModularGCDFiniteField(a, b).multiply(pContentGCD));
// }

for (int uVariable = a.nVariables - 1; uVariable >= 0; --uVariable) {
MultivariatePolynomial<lUnivariatePolynomialZp>
ua = a.asOverUnivariateEliminate(uVariable),
ub = b.asOverUnivariateEliminate(uVariable);

lUnivariatePolynomialZp aContent = ua.content(), bContent = ub.content();
lUnivariatePolynomialZp contentGCD = ua.domain.gcd(aContent, bContent);
lUnivariatePolynomialZp aContent = ua.content(), bContent = ub.content();
lUnivariatePolynomialZp contentGCD = ua.domain.gcd(aContent, bContent);

ua = ua.divideOrNull(aContent);
ub = ub.divideOrNull(bContent);
ua = ua.divideOrNull(aContent);
ub = ub.divideOrNull(bContent);

MultivariatePolynomial<lUnivariatePolynomialZp> ugcd =
ModularGCDFiniteField0(ua, ub, gcdInput.degreeBounds[uVariable], gcdInput.finiteExtensionDegree);
ugcd = ugcd.multiply(contentGCD);
MultivariatePolynomial<lUnivariatePolynomialZp> ugcd =
ModularGCDFiniteField0(ua, ub, gcdInput.degreeBounds[uVariable], gcdInput.finiteExtensionDegree);
if (ugcd == null)
// bad variable choosen
continue;

return gcdInput.restoreGCD(lMultivariatePolynomialZp.asNormalMultivariate(ugcd, uVariable));
ugcd = ugcd.multiply(contentGCD);
return gcdInput.restoreGCD(lMultivariatePolynomialZp.asNormalMultivariate(ugcd, uVariable));
}
throw new RuntimeException();
}

private static final int MAX_OVER_ITERATIONS = 16;
private static final int MAX_OVER_ITERATIONS = 8;

static <uPoly extends IUnivariatePolynomial<uPoly>>
MultivariatePolynomial<uPoly> ModularGCDFiniteField0(
Expand Down Expand Up @@ -1116,7 +1142,15 @@ MultivariatePolynomial<uPoly> ModularGCDFiniteField0(
// over all primes
while (true) {
if (basePrime.degree() >= uDegreeBound + MAX_OVER_ITERATIONS)
throw new RuntimeException();
// fixme:
// probably this should not ever happen, but it happens (extremely rare, only for small
// characteristics and independently on the particular value of MAX_OVER_ITERATIONS)
// the current workaround is to switch to another variable in R[x_N][x1....x_(N-1)]
// representation and try again
//
// UPDATE: when increasing NUMBER_OF_UNDER_DETERMINED_RETRIES the problem seems to be disappeared
// (at the expense of longer time spent in LinZip)
return null;

uPoly prime = primesLoop.next();
fField = new FiniteField<>(prime);
Expand Down Expand Up @@ -1223,10 +1257,12 @@ static <E> MultivariatePolynomial<E> interpolateGCD(MultivariatePolynomial<E> a,
a.checkSameDomainWith(b);
a.checkSameDomainWith(skeleton);

MultivariatePolynomial<E> content = contentGCD(a, b, 0, MultivariateGCD::ZippelGCD);
a = divideExact(a, content);
b = divideExact(b, content);
skeleton = divideSkeletonExact(skeleton, content);
// a and b must be content-free
assert contentGCD(a, b, 0, MultivariateGCD::PolynomialGCD).isConstant();
// lMultivariatePolynomialZp content = contentGCD(a, b, 0, MultivariateGCD::PolynomialGCD);
// a = divideExact(a, content);
// b = divideExact(b, content);
// skeleton = divideSkeletonExact(skeleton, content);

SparseInterpolation<E> interpolation = createInterpolation(-1, a, b, skeleton, rnd);
if (interpolation == null)
Expand All @@ -1235,7 +1271,7 @@ static <E> MultivariatePolynomial<E> interpolateGCD(MultivariatePolynomial<E> a,
if (gcd == null)
return null;

return gcd.multiply(content);
return gcd;//.multiply(content);
}

/* ===================================== Multivariate GCD over finite fields ==================================== */
Expand Down Expand Up @@ -1802,7 +1838,9 @@ public final MultivariatePolynomial<E> evaluate(E newPoint) {
}

/** Number of retries when raise condition occurs; then drop up with new homomorphism */
private static final int NUMBER_OF_UNDER_DETERMINED_RETRIES = 8;
private static final int
NUMBER_OF_UNDER_DETERMINED_RETRIES = 8,
NUMBER_OF_UNDER_DETERMINED_RETRIES_SMALL_FIELD = 24;

static final class LinZipInterpolation<E> extends ASparseInterpolation<E> {
LinZipInterpolation(Domain<E> domain, int variable, MultivariatePolynomial<E> a, MultivariatePolynomial<E> b, Set<DegreeVector> globalSkeleton, TIntObjectHashMap<MultivariatePolynomial<E>> univarSkeleton, int[] sparseUnivarDegrees, int[] evaluationVariables, E[] evaluationPoint, PrecomputedPowersHolder<E> powers, RandomGenerator rnd) {
Expand All @@ -1825,7 +1863,11 @@ public MultivariatePolynomial<E> evaluate() {
int nUnknowns = globalSkeleton.size(), nUnknownScalings = -1;
int raiseFactor = 0;

for (int nTries = 0; nTries < NUMBER_OF_UNDER_DETERMINED_RETRIES; ++nTries) {
// choose dynamically, depending on the cardinality of cyclic group (equal to domain.characteristics)
int nTries = domain.characteristics().bitLength() < 10
? NUMBER_OF_UNDER_DETERMINED_RETRIES_SMALL_FIELD
: NUMBER_OF_UNDER_DETERMINED_RETRIES;
for (int iTry = 0; iTry < nTries; ++iTry) {
int previousFreeVars = -1, underDeterminedTries = 0;
for (; ; ) {
// increment at each loop!
Expand Down Expand Up @@ -1859,7 +1901,7 @@ public MultivariatePolynomial<E> evaluate() {
break;


if (underDeterminedTries > NUMBER_OF_UNDER_DETERMINED_RETRIES) {
if (underDeterminedTries > nTries) {
// raise condition: new equations does not fix enough variables
return null;
}
Expand All @@ -1875,11 +1917,9 @@ public MultivariatePolynomial<E> evaluate() {

MultivariatePolynomial<E> result = a.createZero();
SystemInfo info = solveLinZip(a, systems, nUnknownScalings, result);
if (info == SystemInfo.UnderDetermined) {
if (info == SystemInfo.UnderDetermined)
//try to generate more equations
continue;
}

if (info == SystemInfo.Consistent)
//well done
return result;
Expand Down Expand Up @@ -2744,7 +2784,10 @@ public lMultivariatePolynomialZp evaluate() {
int nUnknowns = globalSkeleton.size(), nUnknownScalings = -1;
int raiseFactor = 0;

for (int nTries = 0; nTries < NUMBER_OF_UNDER_DETERMINED_RETRIES; ++nTries) {
int nTries = domain.modulus < 1024
? NUMBER_OF_UNDER_DETERMINED_RETRIES_SMALL_FIELD
: NUMBER_OF_UNDER_DETERMINED_RETRIES;
for (int iTry = 0; iTry < nTries; ++iTry) {
int previousFreeVars = -1, underDeterminedTries = 0;
for (; ; ) {
// increment at each loop!
Expand Down Expand Up @@ -2778,7 +2821,7 @@ public lMultivariatePolynomialZp evaluate() {
break;


if (underDeterminedTries > NUMBER_OF_UNDER_DETERMINED_RETRIES) {
if (underDeterminedTries > nTries) {
// raise condition: new equations does not fix enough variables
return null;
}
Expand All @@ -2794,11 +2837,9 @@ public lMultivariatePolynomialZp evaluate() {

lMultivariatePolynomialZp result = a.createZero();
SystemInfo info = solveLinZip(a, systems, nUnknownScalings, result);
if (info == SystemInfo.UnderDetermined) {
if (info == SystemInfo.UnderDetermined)
//try to generate more equations
continue;
}

if (info == SystemInfo.Consistent)
//well done
return result;
Expand Down
29 changes: 28 additions & 1 deletion src/test/java/cc/r2/core/poly/multivar/MultivariateGCDTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1270,6 +1270,32 @@ public void testSmallDomain_sparse_variables_random() throws Exception {
GCDAlgorithm.named("Modular gcd (small cardinality)", MultivariateGCD::ModularGCDFiniteField));
}


@Test
public void testSmallDomain2() throws Exception {
lIntegersModulo domain = new lIntegersModulo(3);
String[] vars = {"a", "b", "c", "d", "e"};

lMultivariatePolynomialZp arr[] = {
lMultivariatePolynomialZp.parse("1+2*c^3*d^2+2*b^3*c^3*d^3*e+a*c^3*d*e+2*a^2*b^3*c^2*d^2*e^3+a^2*b^3*c^3*e^2", domain, vars),
lMultivariatePolynomialZp.parse("1+b^3*c^2*d^3*e^3+a*c^3*d*e^2+2*a^3*e^3+2*a^3*b^3*d*e^3+2*a^3*b^3*c*d^3*e", domain, vars),
lMultivariatePolynomialZp.parse("1+2*a*b^3*c+a^2*d^3*e", domain, vars),
lMultivariatePolynomialZp.parse("1+2*b^3*c^3*d^3*e+2*a*b^2*c*d^2*e^3+a*b^3*c^2*d*e^2+a^3*b^2*c^3*d^2", domain, vars),
}, base = arr[0].createOne().multiply(arr);

lMultivariatePolynomialZp a = base;
lMultivariatePolynomialZp b = a.derivative(1);

for (int i = 0; i < its(5, 5); i++) {
timestamp();
lMultivariatePolynomialZp gcd = ModularGCDFiniteField(a, b);
timeElapsed();

assertTrue(dividesQ(a, gcd));
assertTrue(dividesQ(b, gcd));
}
}

@Test
public void testArrayGCD1() throws Exception {
lIntegersModulo domain = new lIntegersModulo(BigPrimes.nextPrime(1321323));
Expand Down Expand Up @@ -2085,7 +2111,8 @@ final GCDSample<Term, Poly> nextSample() {
}

final GCDSample<Term, Poly> nextSample(boolean primitive, boolean monic) {
GCDSample<Term, Poly> sample = nextSample0(primitive, monic);
GCDSample<Term, Poly> sample;
do {sample = nextSample0(primitive, monic);} while (sample.gcd.isZero());
medFactorsDegree.addValue(sample.a.degreeSum());
medFactorsDegree.addValue(sample.b.degreeSum());
medGCDDegree.addValue(sample.gcd.degreeSum());
Expand Down

0 comments on commit 9a0d819

Please sign in to comment.