Skip to content

Commit

Permalink
Major performance improvements for Zippel GCD algorithm for huge poly…
Browse files Browse the repository at this point in the history 
…nomials (#36)

This fixes #35

- efficient methods for evaluation of multivariate polynomials with the use of sparse recursive Horner scheme
- switch between plain and sparse recursive evaluation methods in repeated evaluations in Zippel algorithm
- make some costly things in multivariate polynomials cacheable (first of all degrees())
- use larger primes (around 60 bit) in modular algorithm over Z when the resulting coefficients are expected to be large
  • Loading branch information
@PoslavskySV
PoslavskySV authored Mar 25, 2018
1 parent 219dd3f commit dfeebdf
Show file tree
Hide file tree
Showing 18 changed files with 1,975 additions and 316 deletions.
4 changes: 2 additions & 2 deletions rings.scaladsl/build.sbt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ organization := "cc.redberry"

name := "rings.scaladsl"

version := "2.3"
version := "2.3.1"

scalaVersion := "2.12.3"

Expand All @@ -15,7 +15,7 @@ moduleName := name.value
resolvers += Resolver.mavenLocal

libraryDependencies ++= Seq(
"cc.redberry" % "rings" % "2.3",
"cc.redberry" % "rings" % "2.3.1",
"junit" % "junit" % "4.12" % Test,
"com.novocode" % "junit-interface" % "0.11" % Test exclude("junit", "junit-dep")
)
Expand Down
6 changes: 3 additions & 3 deletions rings/pom.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@

<groupId>cc.redberry</groupId>
<artifactId>rings</artifactId>
<version>2.3</version>
<version>2.3.1</version>
<packaging>jar</packaging>
<name>rings</name>
<url>https://github.com/PoslavskySV/rings/</url>
Expand Down Expand Up @@ -95,12 +95,12 @@
<artifactId>maven-surefire-plugin</artifactId>
<version>2.20.1</version>
<configuration>
<argLine>-Xmx2548m</argLine>
<argLine>-Xmx3060m</argLine>
<trimStackTrace>true</trimStackTrace>
<useFile>false</useFile>
<disableXmlReport>true</disableXmlReport>
<runOrder>random</runOrder>
<!--<forkCount>0</forkCount>-->
<forkCount>1</forkCount>
<reuseForks>false</reuseForks>
</configuration>
</plugin>
Expand Down
156 changes: 100 additions & 56 deletions rings/src/main/java/cc/redberry/rings/poly/multivar/AMultivariatePolynomial.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -290,7 +290,10 @@ public final Poly setOrdering(Comparator<DegreeVector> newOrdering) {
}

/** release caches */
protected void release() { /* add cache in the future */ }
protected void release() {
cachesDegrees = null;
cachedDegree = -1;
}

/**
* Returns the number of terms in this polynomial
Expand Down Expand Up @@ -361,37 +364,6 @@ final Poly loadFrom(MonomialSet<Term> map) {
return self;
}

/**
* Sparsity level: size / (product of degrees)
*/
public double sparsity() {
double sparsity = size();
for (int d : degrees()) {
if (d != 0)
sparsity /= (d + 1);
}
return sparsity;
}

/**
* Sparsity level: {@code size / nDenseTerms} where nDenseTerms is a total number of possible distinct terms with
* total degree not larger than distinct total degrees presented in this.
*/
public double sparsity2() {
TIntHashSet distinctTotalDegrees = new TIntHashSet();
terms.keySet().stream().mapToInt(dv -> dv.totalDegree).forEach(distinctTotalDegrees::add);
TIntIterator it = distinctTotalDegrees.iterator();
double nDenseTerms = 0.0;
while (it.hasNext()) {
int deg = it.next();
double d = BigIntegerUtil.binomial(deg + nVariables - 1, deg).doubleValue();
nDenseTerms += d;
if (d == Double.MAX_VALUE)
return size() / d;
}
return size() / nDenseTerms;
}

/**
* Makes a copy of this with the specified variable dropped
*
Expand Down Expand Up @@ -422,7 +394,7 @@ public final Poly dropVariables(int[] variables) {
*
* @param variables the variables
*/
public final Poly dropSelectVariables(int[] variables) {
public final Poly dropSelectVariables(int... variables) {
MonomialSet<Term> newData = new MonomialSet<>(ordering);
for (Term term : terms)
newData.add(term.dropSelect(variables));
Expand Down Expand Up @@ -488,14 +460,17 @@ final Poly joinNewVariables(int newNVariables, int[] mapping) {
* @return the number of variables (those which are not units)
*/
public final int nUsedVariables() {
int[] degrees = degrees();
int[] degrees = degreesRef();
int r = 0;
for (int d : degrees)
if (d != 0)
++r;
return r;
}

/** cached degree() */
private int cachedDegree = -1;

/**
* Returns the total degree of this polynomial, that is the maximal total degree among all terms
*
Expand All @@ -504,10 +479,13 @@ public final int nUsedVariables() {
@Override
public int degree() {
// fixme replace with degreeSum ?
int max = 0;
for (Term db : terms)
max = Math.max(max, db.totalDegree);
return max;
if (cachedDegree == -1) {
int max = 0;
for (Term db : terms)
max = Math.max(max, db.totalDegree);
cachedDegree = max;
}
return cachedDegree;
}

/**
Expand All @@ -516,7 +494,7 @@ public int degree() {
* @return the maximal degree of variables in this polynomial
*/
public int degreeMax() {
return ArraysUtil.max(degrees());
return ArraysUtil.max(degreesRef());
}

/**
Expand All @@ -526,10 +504,23 @@ public int degreeMax() {
* @return the degree of this polynomial with respect to specified variable
*/
public final int degree(int variable) {
int max = 0;
for (Term db : terms)
max = Math.max(max, db.exponents[variable]);
return max;
return degreesRef()[variable];
}

/** cached degrees */
private int[] cachesDegrees = null;

/** returns reference (content must not be modified) */
protected int[] degreesRef() {
if (cachesDegrees == null) {
int[] degrees = new int[nVariables];
for (Term db : terms)
for (int i = 0; i < nVariables; i++)
if (db.exponents[i] > degrees[i])
degrees[i] = db.exponents[i];
return cachesDegrees = degrees;
}
return cachesDegrees;
}

/**
Expand All @@ -540,12 +531,7 @@ public final int degree(int variable) {
*/
@Override
public final int[] degrees() {
int[] degrees = new int[nVariables];
for (Term db : terms)
for (int i = 0; i < nVariables; i++)
if (db.exponents[i] > degrees[i])
degrees[i] = db.exponents[i];
return degrees;
return degreesRef().clone();
}

/**
Expand Down Expand Up @@ -576,7 +562,7 @@ public final int[] degrees(int variable) {
*/
@Override
public final int degreeSum() {
return ArraysUtil.sum(degrees());
return ArraysUtil.sum(degreesRef());
}

/**
Expand All @@ -586,6 +572,37 @@ public final int totalDegree() {
return degreeSum();
}

/**
* Sparsity level: size / (product of degrees)
*/
public double sparsity() {
double sparsity = size();
for (int d : degreesRef()) {
if (d != 0)
sparsity /= (d + 1);
}
return sparsity;
}

/**
* Sparsity level: {@code size / nDenseTerms} where nDenseTerms is a total number of possible distinct terms with
* total degree not larger than distinct total degrees presented in this.
*/
public double sparsity2() {
TIntHashSet distinctTotalDegrees = new TIntHashSet();
terms.keySet().stream().mapToInt(dv -> dv.totalDegree).forEach(distinctTotalDegrees::add);
TIntIterator it = distinctTotalDegrees.iterator();
double nDenseTerms = 0.0;
while (it.hasNext()) {
int deg = it.next();
double d = BigIntegerUtil.binomial(deg + nVariables - 1, deg).doubleValue();
nDenseTerms += d;
if (d == Double.MAX_VALUE)
return size() / d;
}
return size() / nDenseTerms;
}

/**
* Returns degreeSum - lt().totalDegree
*/
Expand Down Expand Up @@ -641,7 +658,7 @@ public final int univariateVariable() {
return 0;
if (nVariables == 1)
return 0;
int[] degrees = degrees();
int[] degrees = degreesRef();
int var = -1;
for (int i = 0; i < nVariables; i++) {
if (degrees[i] != 0) {
Expand Down Expand Up @@ -779,7 +796,20 @@ Poly asMultivariate(UnivariatePolynomial<Poly> univariate, int uVariable, boolea
* @return multivariate polynomial with coefficients being multivariate polynomials polynomials over {@code
* variables} that is polynomial in R[variables][other_variables]
*/
public abstract MultivariatePolynomial<Poly> asOverMultivariateEliminate(int... variables);
public final MultivariatePolynomial<Poly> asOverMultivariateEliminate(int... variables) {
return asOverMultivariateEliminate(variables, ordering);
}

/**
* Converts this to a multivariate polynomial with coefficients being multivariate polynomials polynomials over
* {@code variables} that is polynomial in R[variables][other_variables]
*
* @param variables the variables to separate
* @param prdering monomial order to use for result
* @return multivariate polynomial with coefficients being multivariate polynomials polynomials over {@code
* variables} that is polynomial in R[variables][other_variables]
*/
public abstract MultivariatePolynomial<Poly> asOverMultivariateEliminate(int[] variables, Comparator<DegreeVector> prdering);

/**
* Convert univariate polynomial over multivariate polynomials to a normal multivariate poly
Expand Down Expand Up @@ -864,6 +894,7 @@ public final Poly multiplyByMonomial(int variable, int exponent) {
for (Term term : oldData)
terms.add(term.set(variable, term.exponents[variable] + exponent));

release();
return self;
}

Expand Down Expand Up @@ -901,6 +932,7 @@ public final Poly setLC(int variable, Poly lc) {
it.remove();
}
terms.putAll(lc.terms);
release();
return self;
}

Expand Down Expand Up @@ -999,7 +1031,7 @@ public final Poly divideDegreeVectorOrNull(DegreeVector monomial) {
/** check whether number of variables is the same */
final void checkSameDomainWith(Term oth) {
if (nVariables != oth.exponents.length)
throw new IllegalArgumentException("Combining multivariate polynomials from different fields");
throw new IllegalArgumentException("Combining multivariate polynomials from different fields: this.nVariables = " + nVariables + " oth.nVariables = " + oth.nVariables());
}

/**
Expand Down Expand Up @@ -1074,6 +1106,17 @@ public final Poly add(Term monomial) {
return self;
}

/**
* Puts {@code monomial} to this polynomial replacing the previous entry if was
*/
public final Poly put(Term monomial) {
checkSameDomainWith(monomial);
terms.add(monomial);
release();
return self;
}


/**
* Subtracts {@code monomial} from this polynomial
*
Expand Down Expand Up @@ -1168,6 +1211,7 @@ public final Poly setAllCoefficientsToUnit() {
Term unit = monomialAlgebra.getUnitTerm(nVariables);
for (Map.Entry<DegreeVector, Term> entry : terms.entrySet())
entry.setValue(unit.setDegreeVector(entry.getKey()));
release();
return self;
}

Expand Down Expand Up @@ -1197,7 +1241,7 @@ public final Set<DegreeVector> getSkeletonExcept(int... variables) {
* @param oth other multivariate polynomial
* @return {@code true} if {@code this} and {@code oth} have the same skeleton and {@code false} otherwise
*/
public final boolean sameSkeletonQ(Poly oth) {
public final boolean sameSkeletonQ(AMultivariatePolynomial oth) {
return getSkeleton().equals(oth.getSkeleton());
}

Expand All @@ -1209,7 +1253,7 @@ public final boolean sameSkeletonQ(Poly oth) {
* @return {@code true} if {@code this} and {@code oth} have the same skeleton with respect to specified {@code
* variables} and {@code false} otherwise
*/
public final boolean sameSkeletonQ(Poly oth, int... variables) {
public final boolean sameSkeletonQ(AMultivariatePolynomial oth, int... variables) {
return getSkeleton(variables).equals(oth.getSkeleton(variables));
}

Expand All @@ -1222,7 +1266,7 @@ public final boolean sameSkeletonQ(Poly oth, int... variables) {
* @return {@code true} if {@code this} and {@code oth} have the same skeleton with respect to all except specified
* {@code variables} and {@code false} otherwise
*/
public final boolean sameSkeletonExceptQ(Poly oth, int... variables) {
public final boolean sameSkeletonExceptQ(AMultivariatePolynomial oth, int... variables) {
return getSkeletonExcept(variables).equals(oth.getSkeletonExcept(variables));
}

Expand Down
15 changes: 4 additions & 11 deletions rings/src/main/java/cc/redberry/rings/poly/multivar/GroebnerBasis.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@
import static cc.redberry.rings.linear.LinearSolver.SystemInfo.*;
import static cc.redberry.rings.poly.multivar.Conversions64bit.*;
import static cc.redberry.rings.poly.multivar.MonomialOrder.GREVLEX;
import static cc.redberry.rings.poly.multivar.MonomialOrder.GRLEX;
import static cc.redberry.rings.poly.multivar.MonomialOrder.isGradedOrder;

/**
* Groebner bases.
Expand Down Expand Up @@ -623,13 +623,6 @@ public int compare(Poly a, Poly b) {
}
}

/** whether monomial order is graded */
static boolean isGradedOrder(Comparator<DegreeVector> monomialOrder) {
return monomialOrder == GREVLEX
|| monomialOrder == GRLEX
|| monomialOrder instanceof GrevLexWithPermutation;
}

/** whether this monomial order is OK for use with homogenization-dehomogenization algorithms */
static boolean isHomogenizationCompatibleOrder(Comparator<DegreeVector> monomialOrder) {
return isGradedOrder(monomialOrder) || monomialOrder == MonomialOrder.LEX;
Expand Down Expand Up @@ -3834,7 +3827,7 @@ List<Poly> solveGB0(List<Poly> generators,
// initial ideal viewed as R[u0, u1, ..., uM][x0, ..., xN]
List<MultivariatePolynomial<Poly>> initialIdeal
= generators
.stream().map(p -> pRing.factory().create(p.terms.values()
.stream().map(p -> pRing.factory().create(p.collection()
.stream().map(t -> new Monomial<>(t, cfRing.factory().create(uTerm.setCoefficientFrom(t))))
.collect(Collectors.toList())))
.collect(Collectors.toList());
Expand Down Expand Up @@ -4005,7 +3998,7 @@ List<Poly> gbSolution(Poly factory, MultivariatePolynomial<Poly>[] gbCandidate)
if (!gbElement.stream().allMatch(Poly::isConstant))
return null;

result.add(factory.create(gbElement.terms.values()
result.add(factory.create(gbElement.collection()
.stream()
.map(t -> t.coefficient.lt().setDegreeVector(t)).collect(Collectors.toList())));
}
Expand Down Expand Up @@ -4098,7 +4091,7 @@ else if (isLinear(eq.reducedEquation)) {

/** whether poly is linear */
static boolean isLinear(AMultivariatePolynomial<? extends AMonomial, ?> poly) {
return poly.terms.values().stream().allMatch(t -> t.totalDegree <= 1);
return poly.collection().stream().allMatch(t -> t.totalDegree <= 1);
}

/** A single equation */
Expand Down
Loading

0 comments on commit dfeebdf

Please sign in to comment.