Skip to content

Instantly share code, notes, and snippets.

@tueda
Last active December 30, 2020 06:00
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save tueda/1066019c2a9823e1450eb979e3793be6 to your computer and use it in GitHub Desktop.
Save tueda/1066019c2a9823e1450eb979e3793be6 to your computer and use it in GitHub Desktop.

Comparison for multivariate polynomial GCD routines

A benchmark result for GCD of random polynomials

  • Desktop PC with i5-3470 @ 3.20GHz, 7.6G RAM
  • Scientific Linux CERN SLC release 6.8 (Carbon)
  • FORM 4.2.0 (Jan 24 2018, v4.2.0-36-ga2acd31) 64-bits
  • Mathematica 11.0.1 for Linux x86 (64-bit)
  • Used the default parameters (--nvars=5, -degree=30, --nterms=50, --coeffpow=50, --nwarmups=10, --nproblems=200)

NOTE: timing distributions can be easily affected by properties of randomly generated polynomials. See this comment.

Display the source blob
Display the rendered blob
Raw
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
*.pyc
__pycache__/
/form/
/argparse.py
*.log
*_vs_*.pdf
*.tmp.*
.deps/
libFermat
all:
python run.py
python plot.py FORM.log Mathematica.log
clean:
rm -rf argparse.py form/ *.pyc __pycache__/ *.log *_vs_*.pdf .deps libFermat
#!/bin/sh
""":" .
exec python "$0" "$@"
"""
import math
import re
import matplotlib; matplotlib.use('Agg') # noqa: E702
import matplotlib.pyplot as plt
import run
argparse = run.argparse
__doc__ = """Plot for multivariate polynomial GCD comparison."""
def parse_log(file):
"""Parse the given log file."""
results = []
for l in file.readlines():
dt, gcd = l.rstrip().split(',')
results.append([float(dt), gcd])
return results
def extract_symbols(expr):
"""Extract symbols from the given text."""
symbols = re.findall('x\d+', expr)
symbols = ['$x_{{{0}}}$'.format(s[1:]) for s in set(symbols)]
return symbols
if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument(
'log1',
type=argparse.FileType('r'),
help='first log file',
)
parser.add_argument(
'log2',
type=argparse.FileType('r'),
help='second log file',
)
opts = parser.parse_args()
name1 = opts.log1.name.replace('.log', '')
name2 = opts.log2.name.replace('.log', '')
res1 = parse_log(opts.log1)
res2 = parse_log(opts.log2)
# Hopefully the following code gives all the symbols in the problem.
symbols = set()
for e in [x for _, x in res1 + res2]:
symbols.update(extract_symbols(e))
symbols = sorted(symbols, key=lambda x: (len(x), x))
xx = [t for t, _ in res1]
yy = [t for t, _ in res2]
tt = [t for t in xx + yy if t != 0]
tmin = 10 ** math.floor(math.log10(min(tt)))
tmax = 10 ** math.ceil(math.log10(max(tt)))
plt.axes().set_aspect('equal')
plt.title('GCD in Z[{0}]'.format(','.join(symbols)))
plt.xlabel('{0} (s)'.format(name1))
plt.ylabel('{0} (s)'.format(name2))
plt.xlim([tmin, tmax])
plt.ylim([tmin, tmax])
plt.grid(True)
plt.xscale('log')
plt.yscale('log')
plt.plot([tmin, tmax], [tmin, tmax], '--', color='gray')
plt.scatter(xx, yy, color='pink', alpha=0.5, edgecolors='red')
plt.savefig('{0}_vs_{1}.pdf'.format(name1, name2))
#!/bin/sh
""":" .
exec python "$0" "$@"
"""
import itertools
import os
import random
import subprocess
import sys
import textwrap
__doc__ = """Comparison for multivariate polynomial GCD routines."""
MAXVARS = 40
random.seed(12345) # fix it for reproducibility
try:
from time import perf_counter
except ImportError:
def perf_counter(): # noqa: D103
import time
return time.time()
def ensure_module(module, package=None):
"""Locate the package if the corresponding module is not found.
Download the package into the current directory if the module is not
available. The package is assumed to be standalone and in pure Python so
that it works just by copying the files.
Parameters
----------
module : str
The name of the required module.
package : str
The name of the package containing the required module.
Returns
-------
mod : module
The loaded module.
Example
-------
>>> ensure_module('form', 'python-form')
"""
try:
imported_module = __import__(module)
except ImportError:
import os
import sys
if package is None:
package = module
sys.stderr.write('Locate {0} package...\n'.format(package))
# Check the PyPI information in JSON.
json_url = 'https://pypi.python.org/pypi/{0}/json'.format(package)
json_file = '{0}.tmp.json'.format(package)
os.system('curl -sL -o {0} {1}'.format(json_file, json_url))
# We need the URL for the wheel.
import json
with open(json_file, 'r') as f:
data = json.load(f)
for u in data['urls']:
if u['packagetype'] == 'bdist_wheel':
wheel_url = u['url']
break
else:
raise ImportError('package ({0}) wheel not found'.format(
package))
# Download the wheel file and expand it.
wheel_file = '{0}.tmp.whl'.format(module)
os.system('curl -sL -o {0} {1}'.format(wheel_file, wheel_url))
os.system('unzip -q {0}'.format(wheel_file, module))
# Delete temporary files.
os.system('rm -fr {0} {1} {2}-*.dist-info'.format(
json_file, wheel_file, package.replace('-', '_')))
# Check the module again.
imported_module = __import__(module)
sys.stderr.write('Done\n')
return imported_module
argparse = ensure_module('argparse')
form = ensure_module('form', 'python-form')
def get_form():
"""Return the singleton for FORM."""
if not hasattr(get_form, '_form'):
f = form.open()
f.write('S x1,...,x{0};'.format(MAXVARS))
f.flush()
get_form._form = f
return get_form._form
def random_poly(nvars, degree, nterms, coeffpow):
"""Return a random polynomial."""
assert 1 <= nvars <= MAXVARS
def random_coeff(m):
while True:
n = random.randint(-10**m, 10**m)
if n != 0:
return n
nterms = max(random.randint(nterms * 2 // 3, nterms), 2)
m = random.randint(1, coeffpow)
flag1 = random.randint(0, 3) == 0 # True for 25%
form = get_form()
while True:
polynomial = ''
for i in range(nterms):
coeff = random_coeff(m)
monomial = ''
sum_degree = degree
if i == 0:
if flag1 and nvars <= sum_degree:
for j in range(nvars):
monomial += '*x{0}'.format(j + 1)
sum_degree -= nvars
xx = list(range(nvars))
random.shuffle(xx)
for j in xx:
p = random.randint(0, sum_degree)
sum_degree -= p
if p != 0:
monomial += '*x{0}^{1}'.format(j + 1, p)
polynomial += '+({0}){1}'.format(coeff, monomial)
form.write('#$a={0};'.format(polynomial))
polynomial = form.read('$a')
if polynomial != '0':
return polynomial
def multiply_poly(p, q):
"""Multiply two polynomials."""
form = get_form()
form.write('#$a={0};#$b={1};#$c=mul_($a,$b);'.format(p, q))
return form.read('$c')
def check_poly(p, q):
"""Check if two polynomials are the same up to an unit."""
form = get_form()
form.write('#$a={0};#$b={1};#$c=$a-$b;#$d=termsin_($c);'.format(p, q))
if form.read('$d') == '0':
return True
form.write('#$c=$a+$b;#$d=termsin_($c);')
if form.read('$d') == '0':
return True
return False
class Problem(object):
"""GCD problem."""
def __init__(self, nvars, degree, nterms, coeffpow):
"""Create a problem."""
self.a = random_poly(nvars, degree, nterms, coeffpow)
self.b = random_poly(nvars, degree, nterms, coeffpow)
self.g = random_poly(nvars, degree, nterms, coeffpow)
self.ag = multiply_poly(self.a, self.g)
self.bg = multiply_poly(self.b, self.g)
class ProblemSet(object):
"""Set of GCD problems."""
def __init__(self, nvars, degree, nterms, coeffpow, nwarmups,
nproblems, logname='problems.log'):
"""Create a set of problems."""
self._nvars = nvars
self._nwarmups = nwarmups
self._nproblems = nproblems
self._problems = [Problem(nvars, degree, nterms, coeffpow)
for _ in range(nwarmups + nproblems)]
with open(logname, 'w') as f:
for i, prob in enumerate(self._problems):
if i >= self._nwarmups:
f.write('gcd_({0},{1})\n'.format(prob.ag, prob.bg))
def __iter__(self):
"""Return the iterator."""
return self._problems.__iter__()
@property
def nvars(self):
"""Return the number of variables."""
return self._nvars
@property
def nwarmups(self):
"""Return the number of warm-up problems."""
return self._nwarmups
@property
def nproblems(self):
"""Return the number of problems."""
return self._nproblems
class Solver(object):
"""GCD solver."""
_name = 'None'
def __iter__(self):
"""Return the iterator."""
return self._result.__iter__()
@property
def name(self):
"""Return the solver name."""
return self._name
def solve(self, problems, logname=None):
"""Solve the given set of problems."""
if logname is None:
logname = self._name + '.log'
self._result = self._solve(problems)
assert len(self._result) == problems.nwarmups + problems.nproblems
with open(logname, 'w') as f:
total_time = 0.0
for i, (dt, gcd) in enumerate(self._result):
if gcd == '1':
sys.stderr.write('warning: gcd = 1\n')
if i >= problems.nwarmups:
f.write('{0},{1}\n'.format(dt, gcd))
total_time += dt
print('{0}: {1:.3f}s'.format(self._name, total_time))
def _solve(self, problems):
results = []
for prob in problems:
t0 = perf_counter()
gcd = self._solve_for(prob.ag, prob.bg)
t1 = perf_counter()
results.append([t1 - t0, gcd])
return results
def _solve_for(self, p, q):
pass
def write_problems(self, filename, problems):
"""Write problems to a file."""
with open(filename, 'w') as f:
for prob in problems:
f.write('{0}\n'.format(prob.ag))
f.write('{0}\n'.format(prob.bg))
def read_answers(self, filename):
"""Read answers from a file."""
with open(filename) as f:
lines = f.readlines()
results = []
for l in lines:
dt, gcd = l.rstrip().split(',')
results.append([float(dt), gcd.replace(' ', '')])
return results
def consistent_with(self, other):
"""Return True if another solver gave the same results."""
gcds1 = [x for _, x in self._result]
gcds2 = [x for _, x in other._result]
for x, y in zip(gcds1, gcds2):
if not check_poly(x, y):
return False
return True
class FormSolver(Solver):
"""GCD by FORM via python-form."""
_name = 'FORM'
def __init__(self):
"""Create a solver."""
self._form = form.open()
self._form.write('S x1,...,x{0};'.format(MAXVARS))
self._form.flush()
def _solve_for(self, p, q):
self._form.write('#$a={0};#$b={1};#$c=gcd_($a,$b);'.format(p, q))
return self._form.read('$c')
class MathematicaSolver(Solver):
"""GCD by Mathematica."""
_name = 'Mathematica'
def _solve(self, problems):
with open('math.tmp.prog', 'w') as f:
f.write('''
n = {0}
'''.format(problems.nwarmups + problems.nproblems))
f.write('''
in = OpenRead["math.tmp.in"]
out = OpenWrite["math.tmp.out"];
Do[
s1 = ReadLine[in];
s2 = ReadLine[in];
p1 = ToExpression[s1];
p2 = ToExpression[s2];
res = Timing[PolynomialGCD[p1, p2]];
WriteString[
out,
ToString[res[[1]] // FortranForm] <> "," <>
ToString[res[[2]] // Expand // InputForm] <> "\n"
];
, n];
Close[in];
Close[out];
Exit[];
''')
self.write_problems('math.tmp.in', problems)
os.system('math <math.tmp.prog >/dev/null 2>&1')
results = self.read_answers('math.tmp.out')
os.system('rm math.tmp.prog math.tmp.in math.tmp.out')
return results
class RingsSolver(Solver):
"""GCD by Rings via rings.repl."""
_name = 'Rings'
def _solve(self, problems):
with open('rings.tmp.prog', 'w') as f:
f.write('''
val N = {0}
implicit val R = MultivariateRing(Z, Array({1}))
'''.format(problems.nwarmups + problems.nproblems,
', '.join('"x{0}"'.format(i + 1)
for i in range(problems.nvars))))
f.write('''
import java.io._
import java.nio.file._
val in = Files.newBufferedReader(Paths.get("rings.tmp.in"))
val out = new PrintWriter(new File("rings.tmp.out"))
for (i <- 1 to N) {
val s1 = in.readLine
val s2 = in.readLine
val p = R(s1)
val q = R(s2)
val t1 = System.nanoTime
val gcd = R.gcd(p, q)
val t2 = System.nanoTime
out.println(((t2 - t1) * 1.0e-9).toString + "," +
R.show(gcd))
}
in.close
out.close
''')
self.write_problems('rings.tmp.in', problems)
os.system('rings.repl <rings.tmp.prog >/dev/null 2>&1')
results = self.read_answers('rings.tmp.out')
os.system('rm rings.tmp.prog rings.tmp.in rings.tmp.out')
return results
class AoutSolver(Solver):
"""GCD by a.out.
The binary ``a.out`` will be executed as
``./a.out <number of variables> <input filename> <output filename>``.
The input file contains an even number of lines. Each of two lines are
text representation of a polynomial of Z[x1,...,xn]. ``a.out`` is expected
to compute the GCD of the two polynomials and write the result to the
output file, where each line must have the format in the following form:
``<timing in seconds>,<GCD of two polynomials>``.
"""
_name = 'a.out'
_binname = 'a.out'
_temp_prefix = 'a.out.tmp'
def _solve(self, problems):
infile = '{0}.in'.format(self._temp_prefix)
outfile = '{0}.out'.format(self._temp_prefix)
self._make()
self.write_problems(infile, problems)
self._run(problems.nvars, infile, outfile)
results = self.read_answers(outfile)
os.system('rm {0}.*'.format(self._temp_prefix))
self._clean()
return results
def _make(self):
pass
def _run(self, nvars, infile, outfile):
# ./a.out <number_of_variables> <input-file> <output-file>
os.system('./{0} {1} {2} {3}'.format(
self._binname,
nvars,
infile,
outfile,
))
def _clean(self):
pass
class FermatSolver(AoutSolver):
"""GCD by Fermat via libFermat."""
_name = 'Fermat'
_binname = 'fer.tmp.prog'
_temp_prefix = 'fer.tmp'
def _make(self):
with open('fer.tmp.prog.cc', 'w') as f:
f.write(textwrap.dedent('''
#include <chrono>
#include <fstream>
#include <iostream>
#include <string>
#include <Fermat.h>
#include <FermatExpression.h>
FermatExpression gcd(const FermatExpression& a, const FermatExpression& b) {
FermatExpression r(a.fer());
(*r.fer())(r.name() + ":=GCD(" + a.name() + "," + b.name() + ")");
return r;
}
int main(int, char* argv[]) {
// a.out <fermat-path> <number_of_variables> <input-file> <output-file>
Fermat fer(argv[1]);
for (int i = 0; i < std::stoi(argv[2]); i++) {
fer.addSymbol("x" + std::to_string(i + 1));
}
std::ifstream infile(argv[3]);
std::ofstream outfile(argv[4]);
std::string line1, line2;
FermatExpression a(&fer), b(&fer), c(&fer);
while (std::getline(infile, line1)) {
std::getline(infile, line2);
a = line1;
b = line2;
auto t0 = std::chrono::system_clock::now();
c = gcd(a, b);
auto t1 = std::chrono::system_clock::now();
std::chrono::duration<double> dt = t1 - t0;
outfile << dt.count() << "," << c.str() << std::endl;
}
}
''')) # noqa
with open('fer.tmp.Makefile', 'w') as f:
f.write(textwrap.dedent('''
CXX = g++
CXXFLAGS = -g -O2 -Wall -Wextra -std=c++11 -pedantic
CPPFLAGS = -IlibFermat/include -IlibFermat/pstreams
DEFS =
LDFLAGS =
DEPDIR = .deps
CXXLINK = $(CXX) $(CXXFLAGS) $(LDFLAGS)
BIN1 = fer.tmp.prog
OBJ1 = fer.tmp.prog.o libFermat/src/Fermat.o libFermat/src/FermatExpression.o
LIB1 =
LINK1 = $(CXXLINK)
BIN = $(BIN1)
OBJ = $(OBJ1)
all: $(BIN)
$(BIN1): $(OBJ1)
\t$(LINK1) -o $(BIN1) $(OBJ1) $(LIB1)
mostlyclean:
\trm -rf $(OBJ)
clean:
\trm -rf $(BIN) $(OBJ) $(DEPDIR) *.gch
.SUFFIXES: .o .cpp
.cpp.o:
\t$(CXX) $(DEFS) $(CPPFLAGS) $(CXXFLAGS) -c -MD -o $@ $< && { [ -d $(DEPDIR) ] || mkdir $(DEPDIR); } && mv $*.d $(DEPDIR)
-include $(DEPDIR)/*.d
''')) # noqa
if os.path.isdir('libFermat'):
os.system('cd libFermat && git pull -q origin')
else:
os.system('git clone -q https://github.com/mprausa/libFermat.git')
os.system('make -s -f fer.tmp.Makefile')
def _run(self, nvars, infile, outfile):
os.system('./fer.tmp.prog {0} {1} {2} {3}'.format(
subprocess.check_output('readlink -f {0}'.format(
subprocess.check_output('command -v fer64', shell=True)
), shell=True).strip(),
nvars,
infile,
outfile,
))
class GiNaCSolver(AoutSolver):
"""GCD by GiNaC."""
_name = 'GiNaC'
_binname = 'ginac.tmp.prog'
_temp_prefix = 'ginac.tmp'
def _make(self):
with open('ginac.tmp.prog.cc', 'w') as f:
f.write(textwrap.dedent('''
#include <chrono>
#include <fstream>
#include <iostream>
#include <string>
#include <ginac/ginac.h>
int main(int, char* argv[]) {
GiNaC::symtab table;
for (int i = 0; i < std::stoi(argv[1]); i++) {
std::string x = "x" + std::to_string(i + 1);
table[x] = GiNaC::symbol(x);
}
GiNaC::parser reader(table);
std::ifstream infile(argv[2]);
std::ofstream outfile(argv[3]);
std::string line1, line2;
while (std::getline(infile, line1)) {
std::getline(infile, line2);
GiNaC::ex a = reader(line1);
GiNaC::ex b = reader(line2);
auto t0 = std::chrono::system_clock::now();
GiNaC::ex c = GiNaC::gcd(a, b);
auto t1 = std::chrono::system_clock::now();
std::chrono::duration<double> dt = t1 - t0;
outfile << dt.count() << "," << c << std::endl;
}
}
''')) # noqa
with open('ginac.tmp.Makefile', 'w') as f:
f.write(textwrap.dedent('''
CXX = g++
CXXFLAGS = -g -O2 -Wall -Wextra -std=c++11 -pedantic
CPPFLAGS =
DEFS =
LDFLAGS =
DEPDIR = .deps
CXXLINK = $(CXX) $(CXXFLAGS) $(LDFLAGS)
BIN1 = ginac.tmp.prog
OBJ1 = ginac.tmp.prog.o
LIB1 = -lginac -lcln
LINK1 = $(CXXLINK)
BIN = $(BIN1)
OBJ = $(OBJ1)
all: $(BIN)
$(BIN1): $(OBJ1)
\t$(LINK1) -o $(BIN1) $(OBJ1) $(LIB1)
mostlyclean:
\trm -rf $(OBJ)
clean:
\trm -rf $(BIN) $(OBJ) $(DEPDIR) *.gch
.SUFFIXES: .o .cpp
.cpp.o:
\t$(CXX) $(DEFS) $(CPPFLAGS) $(CXXFLAGS) -c -MD -o $@ $< && { [ -d $(DEPDIR) ] || mkdir $(DEPDIR); } && mv $*.d $(DEPDIR)
-include $(DEPDIR)/*.d
''')) # noqa
os.system('make -s -f ginac.tmp.Makefile')
if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument(
'--nvars',
default=5,
type=int,
help='set number of variables (default: 5)',
metavar='N',
)
parser.add_argument(
'--degree',
default=30,
type=int,
help='set maximum degree of polynomials (default: 30)',
metavar='N',
)
parser.add_argument(
'--nterms',
default=50,
type=int,
help='set maximum number of terms in a polynomial (default: 50)',
metavar='N',
)
parser.add_argument(
'--coeffpow',
default=50,
type=int,
help='set maximum coefficient as 10^N (default: 50)',
metavar='N',
)
parser.add_argument(
'--nwarmups',
default=10,
type=int,
help='set number of warm-up problems (default: 10)',
metavar='N',
)
parser.add_argument(
'--nproblems',
default=200,
type=int,
help='set number of problems (default: 200)',
metavar='N',
)
parser.add_argument(
'--form',
action='append_const',
const=FormSolver,
help='run FORM',
dest='solvers',
)
parser.add_argument(
'--mathematica',
action='append_const',
const=MathematicaSolver,
help='run Mathematica',
dest='solvers',
)
parser.add_argument(
'--rings',
action='append_const',
const=RingsSolver,
help='run Rings',
dest='solvers',
)
parser.add_argument(
'--fermat',
action='append_const',
const=FermatSolver,
help='run Fermat',
dest='solvers',
)
parser.add_argument(
'--ginac',
action='append_const',
const=GiNaCSolver,
help='run GiNaC',
dest='solvers',
)
parser.add_argument(
'--aout',
action='append_const',
const=AoutSolver,
help='run ./aout',
dest='solvers',
)
opts = parser.parse_args()
problems = ProblemSet(
nvars=opts.nvars,
degree=opts.degree,
nterms=opts.nterms,
coeffpow=opts.coeffpow,
nwarmups=opts.nwarmups,
nproblems=opts.nproblems,
)
solvers = []
if opts.solvers:
for s in opts.solvers:
solvers.append(s())
else:
if os.system('type form >/dev/null 2>&1') == 0:
solvers.append(FormSolver())
if os.system('type math >/dev/null 2>&1') == 0:
solvers.append(MathematicaSolver())
if os.system('type rings.repl >/dev/null 2>&1') == 0:
solvers.append(RingsSolver())
for s in solvers:
s.solve(problems)
failed = False
for s1, s2 in itertools.combinations(solvers, 2):
if not s1.consistent_with(s2):
sys.stderr.write(('Warning: results by {0} and {1} are '
'inconsistent.\n').format(s1.name, s2.name))
failed = True
if failed:
exit(1)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment