"""Functions for generating interesting polynomials, e.g. for benchmarking. """
from sympy.core import Add, Mul, Symbol, sympify, Dummy, symbols
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.core.singleton import S
from sympy.polys.polytools import Poly, PurePoly
from sympy.polys.polyutils import _analyze_gens
from sympy.polys.polyclasses import DMP
from sympy.polys.densebasic import (
dmp_zero, dmp_one, dmp_ground, dmp_normal,
dup_from_raw_dict, dmp_raise, dup_random
)
from sympy.polys.densearith import (
dmp_add_term, dmp_neg, dmp_mul, dmp_sqr
)
from sympy.polys.factortools import (
dup_zz_cyclotomic_poly
)
from sympy.polys.domains import ZZ
from sympy.ntheory import nextprime
from sympy.utilities import cythonized, subsets
@cythonized("n,i")
[docs]def swinnerton_dyer_poly(n, x=None, **args):
"""Generates n-th Swinnerton-Dyer polynomial in `x`. """
if n <= 0:
raise ValueError(
"can't generate Swinnerton-Dyer polynomial of order %s" % n)
if x is not None:
x, cls = sympify(x), Poly
else:
x, cls = Dummy('x'), PurePoly
p, elts = 2, [[x, -sqrt(2)],
[x, sqrt(2)]]
for i in xrange(2, n + 1):
p, _elts = nextprime(p), []
neg_sqrt = -sqrt(p)
pos_sqrt = +sqrt(p)
for elt in elts:
_elts.append(elt + [neg_sqrt])
_elts.append(elt + [pos_sqrt])
elts = _elts
poly = []
for elt in elts:
poly.append(Add(*elt))
if not args.get('polys', False):
return Mul(*poly).expand()
else:
return PurePoly(Mul(*poly), x)
[docs]def cyclotomic_poly(n, x=None, **args):
"""Generates cyclotomic polynomial of order `n` in `x`. """
if n <= 0:
raise ValueError(
"can't generate cyclotomic polynomial of order %s" % n)
poly = DMP(dup_zz_cyclotomic_poly(int(n), ZZ), ZZ)
if x is not None:
poly = Poly.new(poly, x)
else:
poly = PurePoly.new(poly, Dummy('x'))
if not args.get('polys', False):
return poly.as_expr()
else:
return poly
[docs]def symmetric_poly(n, *gens, **args):
"""Generates symmetric polynomial of order `n`. """
gens = _analyze_gens(gens)
if n < 0 or n > len(gens) or not gens:
raise ValueError("can't generate symmetric polynomial of order %s for %s" % (n, gens))
elif not n:
poly = S.One
else:
poly = Add(*[ Mul(*s) for s in subsets(gens, int(n)) ])
if not args.get('polys', False):
return poly
else:
return Poly(poly, *gens)
[docs]def random_poly(x, n, inf, sup, domain=ZZ, polys=False):
"""Return a polynomial of degree ``n`` with coefficients in ``[inf, sup]``. """
poly = Poly(dup_random(n, inf, sup, domain), x, domain=domain)
if not polys:
return poly.as_expr()
else:
return poly
@cythonized("n,i,j")
[docs]def interpolating_poly(n, x, X='x', Y='y'):
"""Construct Lagrange interpolating polynomial for ``n`` data points. """
if isinstance(X, str):
X = symbols("%s:%s" % (X, n))
if isinstance(Y, str):
Y = symbols("%s:%s" % (Y, n))
coeffs = []
for i in xrange(0, n):
numer = []
denom = []
for j in xrange(0, n):
if i == j:
continue
numer.append(x - X[j])
denom.append(X[i] - X[j])
numer = Mul(*numer)
denom = Mul(*denom)
coeffs.append(numer/denom)
return Add(*[ coeff*y for coeff, y in zip(coeffs, Y) ])
@cythonized("n,i")
def fateman_poly_F_1(n):
"""Fateman's GCD benchmark: trivial GCD """
Y = [ Symbol('y_' + str(i)) for i in xrange(0, n + 1) ]
y_0, y_1 = Y[0], Y[1]
u = y_0 + Add(*[ y for y in Y[1:] ])
v = y_0**2 + Add(*[ y**2 for y in Y[1:] ])
F = ((u + 1)*(u + 2)).as_poly(*Y)
G = ((v + 1)*(-3*y_1*y_0**2 + y_1**2 - 1)).as_poly(*Y)
H = Poly(1, *Y)
return F, G, H
@cythonized("n,m,i")
def dmp_fateman_poly_F_1(n, K):
"""Fateman's GCD benchmark: trivial GCD """
u = [K(1), K(0)]
for i in xrange(0, n):
u = [dmp_one(i, K), u]
v = [K(1), K(0), K(0)]
for i in xrange(0, n):
v = [dmp_one(i, K), dmp_zero(i), v]
m = n - 1
U = dmp_add_term(u, dmp_ground(K(1), m), 0, n, K)
V = dmp_add_term(u, dmp_ground(K(2), m), 0, n, K)
f = [[-K(3), K(0)], [], [K(1), K(0), -K(1)]]
W = dmp_add_term(v, dmp_ground(K(1), m), 0, n, K)
Y = dmp_raise(f, m, 1, K)
F = dmp_mul(U, V, n, K)
G = dmp_mul(W, Y, n, K)
H = dmp_one(n, K)
return F, G, H
@cythonized("n,i")
def fateman_poly_F_2(n):
"""Fateman's GCD benchmark: linearly dense quartic inputs """
Y = [ Symbol('y_' + str(i)) for i in xrange(0, n + 1) ]
y_0 = Y[0]
u = Add(*[ y for y in Y[1:] ])
H = Poly((y_0 + u + 1)**2, *Y)
F = Poly((y_0 - u - 2)**2, *Y)
G = Poly((y_0 + u + 2)**2, *Y)
return H*F, H*G, H
@cythonized("n,m,i")
def dmp_fateman_poly_F_2(n, K):
"""Fateman's GCD benchmark: linearly dense quartic inputs """
u = [K(1), K(0)]
for i in xrange(0, n - 1):
u = [dmp_one(i, K), u]
m = n - 1
v = dmp_add_term(u, dmp_ground(K(2), m - 1), 0, n, K)
f = dmp_sqr([dmp_one(m, K), dmp_neg(v, m, K)], n, K)
g = dmp_sqr([dmp_one(m, K), v], n, K)
v = dmp_add_term(u, dmp_one(m - 1, K), 0, n, K)
h = dmp_sqr([dmp_one(m, K), v], n, K)
return dmp_mul(f, h, n, K), dmp_mul(g, h, n, K), h
@cythonized("n,i")
def fateman_poly_F_3(n):
"""Fateman's GCD benchmark: sparse inputs (deg f ~ vars f) """
Y = [ Symbol('y_' + str(i)) for i in xrange(0, n + 1) ]
y_0 = Y[0]
u = Add(*[ y**(n + 1) for y in Y[1:] ])
H = Poly((y_0**(n + 1) + u + 1)**2, *Y)
F = Poly((y_0**(n + 1) - u - 2)**2, *Y)
G = Poly((y_0**(n + 1) + u + 2)**2, *Y)
return H*F, H*G, H
@cythonized("n,i")
def dmp_fateman_poly_F_3(n, K):
"""Fateman's GCD benchmark: sparse inputs (deg f ~ vars f) """
u = dup_from_raw_dict({n + 1: K.one}, K)
for i in xrange(0, n - 1):
u = dmp_add_term([u], dmp_one(i, K), n + 1, i + 1, K)
v = dmp_add_term(u, dmp_ground(K(2), n - 2), 0, n, K)
f = dmp_sqr(
dmp_add_term([dmp_neg(v, n - 1, K)], dmp_one(n - 1, K), n + 1, n, K), n, K)
g = dmp_sqr(dmp_add_term([v], dmp_one(n - 1, K), n + 1, n, K), n, K)
v = dmp_add_term(u, dmp_one(n - 2, K), 0, n - 1, K)
h = dmp_sqr(dmp_add_term([v], dmp_one(n - 1, K), n + 1, n, K), n, K)
return dmp_mul(f, h, n, K), dmp_mul(g, h, n, K), h
# A few useful polynomials from Wang's paper ('78).
f_0 = dmp_normal([
[[1, 2, 3], [2]],
[[3]],
[[4, 5, 6], [1, 2, 1], [1]]
], 2, ZZ)
f_1 = dmp_normal([
[[1, 0], []],
[[1, 0, 1], [20, 30], [1, 10, 0]],
[[1, 0], [30, 20], [1, 10, 1, 610], [20, 230, 300]],
[[1, 10, 0], [30, 320, 200], [600, 6000]]
], 2, ZZ)
f_2 = dmp_normal([
[[1], [1, 0], [1, 0, 0], [1, 0, 0, 0]],
[[]],
[[1], [1, 90], [90, 0]],
[[1, -11], [], [1, -11, 0, 0]],
[[]],
[[1, -11], [90, -990]]
], 2, ZZ)
f_3 = dmp_normal([
[[1], [], []],
[[1, 0, 0, 0, 1]],
[[1, 0], [], [], [1, 0]],
[[1], [1, 0, 0, 0], [], [1, 0, 0, 0, 1, 0], []],
[[1, 0, 0, 0, 1], [1, 0, 0, 0, 1, 1, 0, 0], []],
[[1, 0], [1, 0, 0, 0, 0], []]
], 2, ZZ)
f_4 = dmp_normal([
[[-1, 0], [], [], [], [], [], [], [], []],
[[-1, 0, 0, 0], [], [], [], [], []],
[[-1, 0, 0], [], [], [], [-5], [], [], [], [], [], [], [], []],
[[-1, 0, 0, 0, 0], [], [1, 0, 3, 0], [], [-5, 0, 0], [-1, 0, 0, 0],
[], [], [], []],
[[1, 0, 3, 0, 0, 0], [], [], [-1, 0, 0, 0, 0, 0], []],
[[1, 0, 3, 0, 0], [], [], [-1, 0, 0, 0, 0], [5, 0, 15], [], [], [-
5, 0, 0], [], [], [], []],
[[1, 0, 3, 0, 0, 0, 0], [], [], [-1, 0, 0, 0, 0, 0, 0], [5, 0, 15,
0, 0], [1, 0, 3, 0, 0, 0], [], [-5, 0, 0, 0, 0], []],
[[1, 0, 3, 0, 0, 0, 0, 0]],
[[1, 0, 3, 0, 0, 0, 0], [], [], [], [5, 0, 15, 0, 0], [], [], []],
[[1, 0, 3, 0, 0, 0, 0, 0, 0], [], [], [], [5, 0, 15, 0, 0, 0, 0]]
], 2, ZZ)
f_5 = dmp_normal([
[[-1]],
[[-3], [3, 0]],
[[-3], [6, 0], [-3, 0, 0]],
[[-1], [3, 0], [-3, 0, 0], [1, 0, 0, 0]]
], 2, ZZ)
f_6 = dmp_normal([
[[[2115]], [[]]],
[[[45, 0, 0], [], [], [-45, 0, 0]]],
[[[]]],
[[[-423]], [[-47]], [[]], [[141], [], [94, 0], []], [[]]],
[[[-9, 0, 0], [], [], [9, 0, 0]],
[[-1, 0, 0], [], [], [1, 0, 0]],
[[]],
[[3, 0, 0], [], [2, 0, 0, 0], [-3, 0, 0], [], [-2, 0, 0, 0], []]
]
], 3, ZZ)
w_1 = dmp_normal([
[[4, 0, 0], [4, 0, 0, 0], [-4, 0, 0, 0, 0], [-4, 0, 0, 0, 0, 0], []],
[[1, 0, 0, 0], [12, 0], [-1, 0, 0, 12, 0, 0], [-12, 0, 0, 0], [-12,
0, 0, 0, 0]],
[[8], [6, 8, 0], [-4, 4, -8, 0, 0], [-4, -2, -8, 0, 0, 0], []],
[[2, 0], [1, 0, 0, 0], [-1, 0, -2, 0, 9, 0], [-12, 12, 0, 0], [-12,
3, 0, 0, 0]],
[[6], [-6, 8, 0], [-2, -8, 2, 0, 0], []],
[[2, 0], [-2, 0, 0, 0], [-3, 0], [3, 0, 0, 0]],
[[-2], [2, 0, 0], []]
], 2, ZZ)
w_2 = dmp_normal([
[24, 48, 0, 0],
[24, 0, 0, -72, 0, 0],
[25, 2, 0, 4, 8],
[1, 0, 0, 1, 0, 0, -12],
[1, -1, -2, 292, 0, 0],
[-1, 0, 0, 3, 0, 0, 0],
[-1, 0, 12, 0, 0, 48],
[],
[-12, 0, 0, 0]
], 1, ZZ)