A benchmark for polynomial GCDs

GCD benchmark.md

Tool	Timing (s)
Symbolica	4.1
Mathematica	89
Sympy	3663

mm_gcd.m

a = (1 + 3*x1 + 5*x2 + 7*x3 + 9*x4 + 11*x5 + 13*x6 + 15*x7)^7 - 1
b = (1 - 3*x1 - 5*x2 - 7*x3 + 9*x4 - 11*x5 - 13*x6 + 15*x7)^7 + 1
g = (1 + 3*x1 + 5*x2 + 7*x3 + 9*x4 + 11*x5 + 13*x6 - 15*x7)^7 + 3
ag = Expand[a g]
bg = Expand[b g]
PolynomialGCD[ag, bg] - Expand[g]

symbolica_gcd.py

from symbolica import Expression as E
a = E.parse('(1 + 3*x1 + 5*x2 + 7*x3 + 9*x4 + 11*x5 + 13*x6 + 15*x7)^7 - 1').to_polynomial()
b = E.parse('(1 - 3*x1 - 5*x2 - 7*x3 + 9*x4 - 11*x5 - 13*x6 + 15*x7)^7 + 1').to_polynomial()
g = E.parse('(1 + 3*x1 + 5*x2 + 7*x3 + 9*x4 + 11*x5 + 13*x6 - 15*x7)^7 + 3').to_polynomial()
ag = a * g
bg = b * g
r = ag.gcd(bg) - g

sympy_gcd.py

from sympy.parsing.sympy_parser import parse_expr
from sympy import expand, gcd

a = parse_expr('(1 + 3*x1 + 5*x2 + 7*x3 + 9*x4 + 11*x5 + 13*x6 + 15*x7)**7 - 1')
b = parse_expr('(1 - 3*x1 - 5*x2 - 7*x3 + 9*x4 - 11*x5 - 13*x6 + 15*x7)**7 + 1')
g = parse_expr('(1 + 3*x1 + 5*x2 + 7*x3 + 9*x4 + 11*x5 + 13*x6 - 15*x7)**7 + 3')
ag = a * g
bg = b * g
r = gcd(ag, bg) - expand(g)

Wolfram Mathematica creates ASTs when using expand. That test isn't comparable (they're not apples to apples), among other things because what's the point of using expand in the test? It's to be able to run the test, not for any truly functional purpose.

The timing is solely based on the PolynomialGCD[ag, bg] instruction. The expansions for ag and bg take 1.7 seconds each, and the expansion of g takes only 0.003636s in Mathematica. All three tests work on expanded polynomials, so the comparison is fair, and the polynomials used are the Fateman benchmark polynomials.

Okay, since Mathematica isn't open source, I haven't been able to verify this. Perhaps it would be helpful if the test explained the benchmark for symbolic expansion (with AST tree creation) to clarify things. Otherwise, one might think that the GCD can be fed with polynomials without AST structure versus polynomials with AST structure.

One problem of Mathematica is that you can't really have special structures explicitly, to my knowledge. This is inefficient, but it is also the reality that you have to deal with when you call a function like Cancel on a general expression which uses PolynomialGCD internally. It is not the main bottleneck in this benchmark however. It will be relatively quick to convert an expanded AST to a polynomial, in Symbolica this takes 1.2 seconds for a*g. The main issue is that polynomial arithmetic is slow, because it doesn't use modern algorithms like heap multiplication and likely doesn't use good heuristics.

Perhaps we should ignore and even marginalize all non-open source software these days.