Skip to content

Commit

Permalink
Use EEZ-GCD for very dense GCD problems (#19)
Browse files Browse the repository at this point in the history
- Switch to EEZ-GCD for very dense problems
- Modular GCD in Z with EEZ-GCD for modular images
- Kaltofen&Monagan generic modular algorithm with EEZ-GCD for modular images
  • Loading branch information
PoslavskySV authored Dec 1, 2017
1 parent 55ceccb commit b376a38
Show file tree
Hide file tree
Showing 11 changed files with 854 additions and 128 deletions.
32 changes: 21 additions & 11 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -494,51 +494,61 @@ Index of algorithms implemented in Rings

(the same interpolation techniques as in ZippelGCD is used):

> - [MultivariateGCD.ModularGCD](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/MultivariateGCD.java)
> - [MultivariateGCD.ZippelGCDInZ](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/MultivariateGCD.java)
25. *Kaltofen's & Monagan's generic modular GCD* (\[KalM99\])
25. *Modular GCD over Z (small primes version)*:

> - [MultivariateGCD.ModularGCDInZ](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/MultivariateGCD.java>)
26. *Kaltofen's & Monagan's generic modular GCD with sparse interpolation* (\[KalM99\])

used for computing multivariate GCD over finite fields of very small cardinality:

> - [MultivariateGCD.KaltofenMonaganSparseModularGCDInGF](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/MultivariateGCD.java>)
27. *Kaltofen's & Monagan's generic modular GCD with EEZ-GCD for modular images* (\[KalM99\])

used for computing multivariate GCD over finite fields of very small cardinality:

> - [MultivariateGCD.ModularGCDInGF](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/MultivariateGCD.java)
> - [MultivariateGCD.KaltofenMonaganEEZModularGCDInGF](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/MultivariateGCD.java>)
26. *Multivariate square-free factorization in zero characteristic (Yun's algorithm)* (\[LeeM13\]):
28. *Multivariate square-free factorization in zero characteristic (Yun's algorithm)* (\[LeeM13\]):

> - [MultivariateSquareFreeFactorization.SquareFreeFactorizationYunZeroCharacteristics](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/MultivariateSquareFreeFactorization.java)
27. *Multivariate square-free factorization in non-zero characteristic (Musser's algorithm)* (\[Muss71\], Sec. 8.3 in \[GeCL92\]):
29. *Multivariate square-free factorization in non-zero characteristic (Musser's algorithm)* (\[Muss71\], Sec. 8.3 in \[GeCL92\]):

> - [MultivariateSquareFreeFactorization.SquareFreeFactorizationMusser](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/MultivariateSquareFreeFactorization.java)
> - [MultivariateSquareFreeFactorization.SquareFreeFactorizationMusserZeroCharacteristics](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/MultivariateSquareFreeFactorization.java)
28. *Bernardin's fast dense multivariate Hensel lifting* (\[Bern99\], \[LeeM13\])
30. *Bernardin's fast dense multivariate Hensel lifting* (\[Bern99\], \[LeeM13\])

both for bivariate case (original Bernardin's paper) and multivariate case (Lee thesis) and both with and without precomputed leading coefficients:

> - [multivar.HenselLifting](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/HenselLifting.java)
29. *Sparse Hensel lifting* (\[Kalt85\], \[LeeM13\])
31. *Sparse Hensel lifting* (\[Kalt85\], \[LeeM13\])

> - [multivar.HenselLifting](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/HenselLifting.java)
30. *Fast dense bivariate factorization with recombination* (\[Bern99\], \[LeeM13\]):
32. *Fast dense bivariate factorization with recombination* (\[Bern99\], \[LeeM13\]):

> - [MultivariateFactorization.bivariateDenseFactorSquareFreeInGF](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/MultivariateFactorization.java)
> - [MultivariateFactorization.bivariateDenseFactorSquareFreeInZ](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/MultivariateFactorization.java)
31. *Kaltofen's multivariate factorization in finite fields* (\[Kalt85\], \[LeeM13\])
33. *Kaltofen's multivariate factorization in finite fields* (\[Kalt85\], \[LeeM13\])

modified version of original Kaltofen's algorithm for leading coefficient precomputation with square-free decomposition (instead of distinct variables decomposition) due to Lee is used; further adaptations are made to work in finite fields of very small cardinality; the resulting algorithm is close to \[LeeM13\], but at the same time has many differences in details:

> - [MultivariateFactorization.factorInGF](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/MultivariateFactorization.java)
32. *Kaltofen's multivariate factorization Z* (\[Kalt85\], \[LeeM13\])
34. *Kaltofen's multivariate factorization Z* (\[Kalt85\], \[LeeM13\])

(with the same modifications as for algorithm for finite fields):

> - [MultivariateFactorization.factorInZ](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/MultivariateFactorization.java)
33. *Multivariate polynomial interpolation with Newton method*:
35. *Multivariate polynomial interpolation with Newton method*:

> - [MultivariateInterpolation](https://github.com/PoslavskySV/rings/tree/develop/rings/src/main/java/cc/redberry/rings/poly/multivar/MultivariateInterpolation.java)
Expand Down
10 changes: 6 additions & 4 deletions doc/guide.rst
Original file line number Diff line number Diff line change
Expand Up @@ -1796,12 +1796,14 @@ Multivariate GCD

- ``BrownGCD`` |br| Brown's GCD for multivariate polynomials over finite fields (see [Brow71]_, Sec 7.4 in [GeCL92]_, [Yang09]_).
- ``ZippelGCD`` |br| Zippel's sparse algorithm for multivariate GCD over fields. Works both in case of monic polynomials with fast Vandermonde linear systems (see [Zipp79]_, [Zipp93]_) and in case of non-monic input (LINZIP, see [dKMW05]_, [Yang09]_).
- ``ModularGCD`` |br| Modular GCD for multivariate polynomials over Z with sparse interpolation (see [Zipp79]_, [Zipp93]_, [dKMW05]_) (the same interpolation techniques as in ZippelGCD is used).
- ``ModularGCDInGF`` |br| Kaltofen's & Monagan's generic modular GCD (see [KalM99]_) for multivariate polynomials over finite fields with very small cardinality.
- ``ZippelGCDInZ`` |br| Zippel's sparse algorithm for multivariate GCD over Z (see [Zipp79]_, [Zipp93]_, [dKMW05]_) (the same interpolation techniques as in ZippelGCD is used)..
- ``ModularGCDInZ`` |br| Standard modular algorithm (small primes) for GCD over Z.
- ``KaltofenMonaganSparseModularGCDInGF`` |br| Kaltofen's & Monagan's generic modular GCD (see [KalM99]_) for multivariate polynomials over finite fields with very small cardinality with sparse Zippel's techniques similar to ZippelGCDInZ
- ``KaltofenMonaganEEZModularGCDInGF`` |br| Kaltofen's & Monagan's generic modular GCD (see [KalM99]_) for multivariate polynomials over finite fields with very small cardinality with EEZ-GCD used for modular images
- ``EZGCD`` |br| Extended Zassenhaus GCD (EZ-GCD) for multivariate polynomials over finite fields (see Sec. 7.6 in [GeCL92]_ and [MosY73]_).
- ``EEZGCD`` |br| Enhanced Extended Zassenhaus GCD (EEZ-GCD) for multivariate polynomials over finite fields (see [Wang80]_).

The upper-level method ``MultivariateGCD.PolynomialGCD`` uses ``ZippelGCD`` for polynomials over finite fields (it shows the best performance in practice). In case of finite fields of very small cardinality ``ModularGCDInGF`` is used. ``ModularGCD`` is used for polynomials in :math:`Z[x]` and :math:`Q[x]`. Algorithms ``BrownGCD`` and ``ZippelGCD`` automatically switch to ``ModularGCDInGF`` in case if the coefficient domain has insufficiently large cardinality. Examples:
The upper-level method ``MultivariateGCD.PolynomialGCD`` switches between Zippel-like algorithms and EEZ-GCD based algorithms. The latter are used only on a very dense problems (which occur rarely), while the former are actually used in most cases. In case of finite fields of very small cardinality Kaltofen's & Monagan's algorithm is used. Examples:

.. tabs::

Expand Down Expand Up @@ -1950,7 +1952,7 @@ Details of implementation can be found in `MultivariateGCD`_.
.. _ref-multivariate-factorization:

Multivariate factorization
^^^^^^^^^^^^^^^^^^^^^^^^^^
^^^^^^^^^^^^^^^^^^^^^^^^^^

Implementation of multivariate factorization in |Rings| is distributed over two classes:

Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,84 @@
package cc.redberry.rings.benchmark;

import cc.redberry.rings.Ring;
import cc.redberry.rings.Rings;
import cc.redberry.rings.bigint.BigInteger;
import cc.redberry.rings.poly.PolynomialMethods;
import cc.redberry.rings.poly.multivar.AMultivariatePolynomial;
import cc.redberry.rings.poly.multivar.MultivariateGCD;
import cc.redberry.rings.poly.multivar.MultivariatePolynomial;
import cc.redberry.rings.util.TimeUnits;

import java.util.Arrays;

/**
* @author Stanislav Poslavsky
* @since 2.2
*/
public class Rings_vs_Singular_vs_Mathematica_GCD2 {
public static void main(String[] args) throws Exception {

for (Ring<BigInteger> ring : new Ring[]{Rings.Z, Rings.Zp((1 << 19) - 1)}) {
System.out.println("Ring: " + ring);

for (int exp = 4; exp < 7; exp++) {
System.out.println("\n\n\n\n");
System.out.println("<><><><><><><><><><><><><><><><><><><><><><><><><><><>");
System.out.println("exp: " + exp);

MultivariatePolynomial<BigInteger>
a = MultivariatePolynomial.parse("1 + 3*a + 5*b + 7*c + 9*d + 11*e + 13*f + 15*g", ring),
b = MultivariatePolynomial.parse("1 - 3*a + 5*b - 7*c + 9*d - 11*e + 13*f - 15*g", ring),
g = MultivariatePolynomial.parse("1 + 3*a - 5*b + 7*c - 9*d + 11*e - 13*f + 15*g", ring);

a = PolynomialMethods.polyPow(a, exp);
a.decrement();
b = PolynomialMethods.polyPow(b, exp);
g = PolynomialMethods.polyPow(g, exp);
g.increment();

MultivariatePolynomial<BigInteger> ag = a.clone().multiply(g);
MultivariatePolynomial<BigInteger> bg = b.clone().multiply(g);

// info(a, "a");
// info(b, "b");
// info(g, "g");
//
// info(ag, "ag");
// info(bg, "bg");

System.out.println("\n=================\n");
for (int i = 0; i < 5; ++i) {
long start = System.nanoTime();
int size = MultivariateGCD.PolynomialGCD(ag, bg).size();
long ringsTime = System.nanoTime() - start;

System.out.println("Rings: " + TimeUnits.nanosecondsToString(ringsTime));

long singularTime = Bench.singularGCD(ag, bg).nanoTime;
System.out.println("Singular: " + TimeUnits.nanosecondsToString(singularTime));

long mmaTime = Bench.mathematicaGCD(ag, bg).nanoTime;
System.out.println("MMA: " + TimeUnits.nanosecondsToString(mmaTime));

System.out.println();
}
}
}
}

private static void info(AMultivariatePolynomial p) {
info(p, "");
}

private static void info(AMultivariatePolynomial p, String poly) {
System.out.println();
if (!poly.isEmpty())
System.out.println("Info for " + poly);
System.out.println("Ring: " + p.coefficientRingToString());
System.out.println("Degrees: " + Arrays.toString(p.degrees()));
System.out.println("Size: " + p.size());
System.out.println("Sparsity: " + p.sparsity());
System.out.println("Sparsity2: " + p.sparsity2());
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -214,10 +214,12 @@ abstract class PolynomialRing[Poly <: IPolynomial[Poly], E]
protected[scaladsl] final def divRem(a: Poly, b: Poly): (Poly, Poly) = {
val qd = theRing.divideAndRemainder(a, b)
if (qd == null)
throw new ArithmeticException(s"not divisible with remainder: ${this show a} / ${this show b}")
throw new ArithmeticException(s"not divisible with remainder: ${this show a } / ${this show b }")
(qd(0), qd(1))
}

final def apply(value: E): Poly = getConstant(value)

/**
* Constant polynomial with specified value
*/
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -795,4 +795,102 @@ class Examples {
println(xan)
println(xna)
}

@Test
def test30: Unit = {
import syntax._

/**
* Some generic function
*
* @tparam Poly type of polynomial
* @tparam Coef type of polynomial coefficients
*/
def genericFunc[Poly <: IPolynomial[Poly], Coef]
(poly: Poly)(implicit ring: PolynomialRing[Poly, Coef]): Poly = {
// implicit coefficient ring (setups algebraic operators on type Coef)
implicit val cfRing: Ring[Coef] = ring.cfRing

// c.c. and l.c. are of type Coef
val (cc: Coef, lc: Coef) = (poly.cc, poly.lc)
// do some calculations
val l: Coef = (cc + lc).pow(10)
// return the result
if (ring.characteristic === 2)
poly - l
else
l + l * poly
}

implicit val uRing = UnivariateRing(Z, "x")
// Coef will be inferred as IntZ and Poly as UnivariatePolynomial[IntZ]
val ur: UnivariatePolynomial[IntZ] = genericFunc(uRing("x^2 + 2"))

implicit val mRing = MultivariateRingZp64(17, Array("x", "y", "z"))
// Coef will be inferred as Long and Poly as MultivariatePolynomialZp64
val mr: MultivariatePolynomialZp64 = genericFunc(mRing("x^2 + 2 + y + z"))
}


@Test
def test31: Unit = {
import syntax._

/** Lagrange polynomial interpolation formula */
def lagrange[Poly <: IUnivariatePolynomial[Poly], Coef]
(points: Seq[Coef], values: Seq[Coef])
(implicit ring: IUnivariateRing[Poly, Coef]) = {
// implicit coefficient ring (setups algebraic operators on type Coef)
implicit val cfRing: Ring[Coef] = ring.cfRing
if (!cfRing.isField) throw new IllegalArgumentException
points.indices
.foldLeft(ring(0)) { case (sum, i) =>
sum + points.indices
.filter(_ != i)
.foldLeft(ring(values(i))) { case (product, j) =>
product * (ring.`x` - points(j)) / (points(i) - points(j))
}
}
}

def interpolate[Poly <: IUnivariatePolynomial[Poly], Coef]
(points: Seq[(Coef, Coef)])
(implicit ring: IUnivariateRing[Poly, Coef]) = {
// implicit coefficient ring (setups algebraic operators on type Coef)
implicit val cfRing: Ring[Coef] = ring.cfRing
if (!cfRing.isField) throw new IllegalArgumentException
points.indices
.foldLeft(ring(0)) { case (sum, i) =>
sum + points.indices
.filter(_ != i)
.foldLeft(ring(points(i)._2)) { case (product, j) =>
product * (ring.`x` - points(j)._1) / (points(i)._1 - points(j)._1)
}
}
}

// coefficient ring Frac(Z/13[a,b,c])
val cfRing = Frac(MultivariateRingZp64(2, Array("a", "b", "c")))
val (a, b, c) = cfRing("a", "b", "c")

implicit val ring = UnivariateRing(cfRing, "x")
// interpolate with Lagrange formula
val data = Seq(a -> b, b -> c, c -> a)
val poly = interpolate(data)
assert(data.forall { case (p, v) => poly.eval(p) == v })

println(poly)


//UnivariatePolynomial[Rational[MultivariatePolynomialZp64]]

// // Q[x]
// implicit val ringA = UnivariateRing(Q, "x")
// // Coef will be inferred as IntZ and Poly as UnivariatePolynomial[IntZ]
// val ur: UnivariatePolynomial[Rational[IntZ]]
// = lagrange(
// Array(1, 2, 3).map(Q(_)),
// Array(5, 4, 3).map(Q(_))
// )
}
}
24 changes: 23 additions & 1 deletion rings/src/main/java/cc/redberry/rings/bigint/BigIntegerUtil.java
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,7 @@ public static BigInteger pow(final BigInteger base, long exponent) {

if (exponent < Integer.MAX_VALUE)
return pow(base, (int) exponent);

BigInteger result = ONE;
BigInteger k2p = base;
for (; ; ) {
Expand Down Expand Up @@ -210,4 +210,26 @@ else if (ab.compareTo(n) < 0)
}
return null;
}

/**
* Factorial of a number
*/
public static BigInteger factorial(int number) {
BigInteger r = BigInteger.ONE;
for (int i = 1; i <= number; i++)
r = r.multiply(i);
return r;
}

/**
* Binomial coefficient
*/
public static BigInteger binomial(int n, int k) {
if (k > n - k)
k = n - k;
BigInteger b = BigInteger.ONE;
for (int i = 1, m = n; i <= k; i++, m--)
b = b.multiply(m).divideExact(BigInteger.valueOf(i));
return b;
}
}
Original file line number Diff line number Diff line change
@@ -1,12 +1,14 @@
package cc.redberry.rings.poly.multivar;

import cc.redberry.rings.WithVariables;
import cc.redberry.rings.bigint.BigIntegerUtil;
import cc.redberry.rings.poly.IPolynomial;
import cc.redberry.rings.poly.MultivariateRing;
import cc.redberry.rings.poly.univar.IUnivariatePolynomial;
import cc.redberry.rings.poly.univar.UnivariatePolynomial;
import cc.redberry.rings.poly.univar.UnivariatePolynomialZp64;
import cc.redberry.rings.util.ArraysUtil;
import gnu.trove.iterator.TIntIterator;
import gnu.trove.set.hash.TIntHashSet;
import org.apache.commons.math3.random.RandomGenerator;

Expand Down Expand Up @@ -314,6 +316,25 @@ public double sparsity() {
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
Loading

0 comments on commit b376a38

Please sign in to comment.