Permalink
Browse files

Cythonized zero-level algorithms in new polynomials module

Pure Python:

In [1]: f = expand(((x+y+z)**15+1)*((x+y+z)**15+2))

In [2]: %time a = factor(f)
CPU times: user 109.45 s, sys: 0.01 s, total: 109.47 s
Wall time: 110.83 s

In [4]: %time a = factor(f)
CPU times: user 109.31 s, sys: 0.03 s, total: 109.34 s
Wall time: 110.68 s

Pure mode Cython:

In [1]: f = expand(((x+y+z)**15+1)*((x+y+z)**15+2))

In [2]: %time a = factor(f)
CPU times: user 72.09 s, sys: 1.02 s, total: 73.11 s
Wall time: 74.18 s

In [4]: %time a = factor(f)
CPU times: user 72.81 s, sys: 0.04 s, total: 72.85 s
Wall time: 73.74 s

On average Cython version is two times faster than pure
Python. This is an improvement and hopefully it should
get even better in future.

To make cooperation with Cython more comfortable a new
decorator was added to sympy/utilities/cythonutils.py.

Example:

@cythonized('n,k')
def some_function(n):
    return n**k

This means that `some_function` will be compiled by
Cython if available, treating variables n and k as
integers (cython.int). This is convenient because
currently we don't use any other types.

If Cython is not available then @cythonized is an
empty decorator (there is no performance penalty).

To take advantage of pure mode Cython, you have to
compile modules which support cythonization. To do
this, issue:

python build.py build_ext --inplace

in SymPy's root directory (or use make). Then run
isympy or import sympy as usually (compiled modules
will have priority over pure Python).
  • Loading branch information...
1 parent 1c43497 commit afb38863484d13b61f0d636bec980f1c0c1fed24 @mattpap mattpap committed Aug 19, 2009
View
9 Makefile
@@ -0,0 +1,9 @@
+# build or clean Cythonized modules
+
+all:
+ python build.py build_ext --inplace
+
+clean:
+ rm -rf sympy/polys/*.{pyc,c,so}
+ rm -rf build/
+
View
50 build.py
@@ -0,0 +1,50 @@
+#!/usr/bin/env python
+
+import os
+
+from Cython.Compiler.Main import compile
+
+from distutils.core import setup, Extension
+from distutils.command.build_ext import build_ext
+
+source_root = os.path.dirname(__file__)
+
+compiled_modules = [
+ "sympy.polys.densearith",
+ "sympy.polys.densebasic",
+ "sympy.polys.densetools",
+ "sympy.polys.galoistools",
+ "sympy.polys.factortools",
+ "sympy.polys.specialpolys",
+ "sympy.polys.monomialtools",
+]
+
+extensions = []
+
+for module in compiled_modules:
+ source_file = os.path.join(source_root, *module.split('.')) + ".py"
+
+ print("Compiling module %s ..." % module)
+ result = compile(source_file)
+
+ if result.c_file is None:
+ raise RuntimeError("failed to compile %s" % module)
+
+ extensions.append(
+ Extension(module, sources=[result.c_file],
+ extra_compile_args=['-O2', '-Wall'],
+ )
+ )
+
+setup(
+ name = "SymPy",
+ packages = [
+ "sympy",
+ "sympy.polys",
+ ],
+ cmdclass = {
+ "build_ext": build_ext
+ },
+ ext_modules = extensions
+)
+
View
44 sympy/polys/densearith.py
@@ -6,13 +6,16 @@
dup_strip, dmp_strip,
dmp_zero_p, dmp_zero,
dmp_one_p, dmp_one,
- dmp_zeros,
+ dmp_zeros
)
from sympy.polys.polyerrors import (
- ExactQuotientFailed,
+ ExactQuotientFailed
)
+from sympy.utilities import cythonized
+
+@cythonized("i,n,m")
def dup_add_term(f, c, i, K):
"""Add `c*x**i` to `f` in `K[x]`. """
if not c:
@@ -29,6 +32,7 @@ def dup_add_term(f, c, i, K):
else:
return f[:m] + [f[m]+c] + f[m+1:]
+@cythonized("i,u,v,n,m")
def dmp_add_term(f, c, i, u, K):
"""Add `c(x_2..x_u)*x_0**i` to `f` in `K[X]`. """
if not u:
@@ -50,6 +54,7 @@ def dmp_add_term(f, c, i, u, K):
else:
return f[:m] + [dmp_add(f[m], c, v, K)] + f[m+1:]
+@cythonized("i,n,m")
def dup_sub_term(f, c, i, K):
"""Subtract `c*x**i` from `f` in `K[x]`. """
if not c:
@@ -66,6 +71,7 @@ def dup_sub_term(f, c, i, K):
else:
return f[:m] + [f[m]-c] + f[m+1:]
+@cythonized("i,u,v,n,m")
def dmp_sub_term(f, c, i, u, K):
"""Subtract `c(x_2..x_u)*x_0**i` from `f` in `K[X]`. """
if not u:
@@ -87,13 +93,15 @@ def dmp_sub_term(f, c, i, u, K):
else:
return f[:m] + [dmp_sub(f[m], c, v, K)] + f[m+1:]
+@cythonized("i")
def dup_mul_term(f, c, i, K):
"""Multiply `f` by `c*x**i` in `K[x]`. """
if not c or not f:
return []
else:
return [ cf * c for cf in f ] + [K.zero]*i
+@cythonized("i,u,v")
def dmp_mul_term(f, c, i, u, K):
"""Multiply `f` by `c(x_2..x_u)*x_0**i` in `K[X]`. """
if not u:
@@ -115,6 +123,7 @@ def dup_mul_ground(f, c, K):
else:
return [ cf * c for cf in f ]
+@cythonized("u,v")
def dmp_mul_ground(f, c, u, K):
"""Multiply `f` by a constant value in `K[X]`. """
if not u:
@@ -133,6 +142,7 @@ def dup_quo_ground(f, c, K):
return [ K.quo(cf, c) for cf in f ]
+@cythonized("u,v")
def dmp_quo_ground(f, c, u, K):
"""Quotient by a constant in `K[X]`. """
if not u:
@@ -154,6 +164,7 @@ def dup_exquo_ground(f, c, K):
else:
return [ cf // c for cf in f ]
+@cythonized("u,v")
def dmp_exquo_ground(f, c, u, K):
"""Exact quotient by a constant in `K[X]`. """
if not u:
@@ -167,6 +178,7 @@ def dup_abs(f, K):
"""Make all coefficients positive in `K[x]`. """
return [ K.abs(coeff) for coeff in f ]
+@cythonized("u,v")
def dmp_abs(f, u, K):
"""Make all coefficients positive in `K[X]`. """
if not u:
@@ -180,6 +192,7 @@ def dup_neg(f, K):
"""Negate a polynomial in `K[x]`. """
return [ -coeff for coeff in f ]
+@cythonized("u,v")
def dmp_neg(f, u, K):
"""Negate a polynomial in `K[X]`. """
if not u:
@@ -189,6 +202,7 @@ def dmp_neg(f, u, K):
return [ dmp_neg(cf, u-1, K) for cf in f ]
+@cythonized("df,dg,k")
def dup_add(f, g, K):
"""Add dense polynomials in `K[x]`. """
if not f:
@@ -211,6 +225,7 @@ def dup_add(f, g, K):
return h + [ a + b for a, b in zip(f, g) ]
+@cythonized("u,v,df,dg,k")
def dmp_add(f, g, u, K):
"""Add dense polynomials in `K[X]`. """
if not u:
@@ -240,6 +255,7 @@ def dmp_add(f, g, u, K):
return h + [ dmp_add(a, b, v, K) for a, b in zip(f, g) ]
+@cythonized("df,dg,k")
def dup_sub(f, g, K):
"""Subtract dense polynomials in `K[x]`. """
if not f:
@@ -262,6 +278,7 @@ def dup_sub(f, g, K):
return h + [ a - b for a, b in zip(f, g) ]
+@cythonized("u,v,df,dg,k")
def dmp_sub(f, g, u, K):
"""Subtract dense polynomials in `K[X]`. """
if not u:
@@ -295,6 +312,7 @@ def dup_add_mul(f, g, h, K):
"""Returns `f + g*h` where `f, g, h` are in `K[x]`. """
return dup_add(f, dup_mul(g, h, K), K)
+@cythonized("u")
def dmp_add_mul(f, g, h, u, K):
"""Returns `f + g*h` where `f, g, h` are in `K[X]`. """
return dmp_add(f, dmp_mul(g, h, u, K), u, K)
@@ -303,10 +321,12 @@ def dup_sub_mul(f, g, h, K):
"""Returns `f - g*h` where `f, g, h` are in `K[x]`. """
return dup_sub(f, dup_mul(g, h, K), K)
+@cythonized("u")
def dmp_sub_mul(f, g, h, u, K):
"""Returns `f - g*h` where `f, g, h` are in `K[X]`. """
return dmp_sub(f, dmp_mul(g, h, u, K), u, K)
+@cythonized("df,dg,i,j")
def dup_mul(f, g, K):
"""Multiply dense polynomials in `K[x]`. """
if f == g:
@@ -330,6 +350,7 @@ def dup_mul(f, g, K):
return h
+@cythonized("u,v,df,dg,i,j")
def dmp_mul(f, g, u, K):
"""Multiply dense polynomials in `K[X]`. """
if not u:
@@ -360,6 +381,7 @@ def dmp_mul(f, g, u, K):
return h
+@cythonized("df,jmin,jmax,n,i,j")
def dup_sqr(f, K):
"""Square dense polynomials in `K[x]`. """
df, h = dup_degree(f), []
@@ -387,6 +409,7 @@ def dup_sqr(f, K):
return h
+@cythonized("u,v,df,jmin,jmax,n,i,j")
def dmp_sqr(f, u, K):
"""Square dense polynomials in `K[X]`. """
if not u:
@@ -422,6 +445,7 @@ def dmp_sqr(f, u, K):
return h
+@cythonized("n,m")
def dup_pow(f, n, K):
"""Raise f to the n-th power in `K[x]`. """
if not n:
@@ -446,6 +470,7 @@ def dup_pow(f, n, K):
return g
+@cythonized("u,n,m")
def dmp_pow(f, n, u, K):
"""Raise f to the n-th power in `K[X]`. """
if not u:
@@ -473,6 +498,7 @@ def dmp_pow(f, n, u, K):
return g
+@cythonized("df,dg,dr,N,j")
def dup_pdiv(f, g, K):
"""Polynomial pseudo-division in `K[x]`. """
df = dup_degree(f)
@@ -511,6 +537,7 @@ def dup_pdiv(f, g, K):
return q, r
+@cythonized("df,dg,dr,N,j")
def dup_prem(f, g, K):
"""Polynomial pseudo-remainder in `K[x]`. """
df = dup_degree(f)
@@ -554,6 +581,7 @@ def dup_pexquo(f, g, K):
"""Polynomial exact pseudo-quotient in `K[X]`. """
return dup_pdiv(f, g, K)[0]
+@cythonized("u,df,dg,dr,N,j")
def dmp_pdiv(f, g, u, K):
"""Polynomial pseudo-division in `K[X]`. """
if not u:
@@ -596,6 +624,7 @@ def dmp_pdiv(f, g, u, K):
return q, r
+@cythonized("u,df,dg,dr,N,j")
def dmp_prem(f, g, u, K):
"""Polynomial pseudo-remainder in `K[X]`. """
if not u:
@@ -645,6 +674,7 @@ def dmp_pexquo(f, g, u, K):
"""Polynomial exact pseudo-quotient in `K[X]`. """
return dmp_pdiv(f, g, u, K)[0]
+@cythonized("df,dg,dr,j")
def dup_rr_div(f, g, K):
"""Univariate division with remainder over a ring. """
df = dup_degree(f)
@@ -680,6 +710,7 @@ def dup_rr_div(f, g, K):
return q, r
+@cythonized("u,df,dg,dr,j")
def dmp_rr_div(f, g, u, K):
"""Multivariate division with remainder over a ring. """
if not u:
@@ -720,6 +751,7 @@ def dmp_rr_div(f, g, u, K):
return q, r
+@cythonized("df,dg,dr,j")
def dup_ff_div(f, g, K):
"""Polynomial division with remainder over a field. """
df = dup_degree(f)
@@ -752,6 +784,7 @@ def dup_ff_div(f, g, K):
return q, r
+@cythonized("u,df,dg,dr,j")
def dmp_ff_div(f, g, u, K):
"""Polynomial division with remainder over a field. """
if not u:
@@ -816,17 +849,20 @@ def dup_exquo(f, g, K):
"""Returns exact polynomial quotient in `K[x]`. """
return dup_div(f, g, K)[0]
+@cythonized("u")
def dmp_div(f, g, u, K):
"""Polynomial division with remainder in `K[X]`. """
if K.has_Field:
return dmp_ff_div(f, g, u, K)
else:
return dmp_rr_div(f, g, u, K)
+@cythonized("u")
def dmp_rem(f, g, u, K):
"""Returns polynomial remainder in `K[X]`. """
return dmp_div(f, g, u, K)[1]
+@cythonized("u")
def dmp_quo(f, g, u, K):
"""Returns polynomial quotient in `K[X]`. """
q, r = dmp_div(f, g, u, K)
@@ -836,6 +872,7 @@ def dmp_quo(f, g, u, K):
else:
raise ExactQuotientFailed('%s does not divide %s' % (g, f))
+@cythonized("u")
def dmp_exquo(f, g, u, K):
"""Returns exact polynomial quotient in `K[X]`. """
return dmp_div(f, g, u, K)[0]
@@ -847,6 +884,7 @@ def dup_max_norm(f, K):
else:
return max(dup_abs(f, K))
+@cythonized("u,v")
def dmp_max_norm(f, u, K):
"""Returns maximum norm of a polynomial in `K[X]`. """
if not u:
@@ -863,6 +901,7 @@ def dup_l1_norm(f, K):
else:
return sum(dup_abs(f, K))
+@cythonized("u,v")
def dmp_l1_norm(f, u, K):
"""Returns l1 norm of a polynomial in `K[X]`. """
if not u:
@@ -884,6 +923,7 @@ def dup_expand(polys, K):
return f
+@cythonized("u")
def dmp_expand(polys, u, K):
"""Multiply together several polynomials in `K[X]`. """
if not polys:
View
355 sympy/polys/densebasic.py
@@ -4,9 +4,11 @@
from sympy.core import igcd, ilcm
from sympy.polys.monomialtools import (
- monomial_min, monomial_div,
+ monomial_min, monomial_div
)
+from sympy.utilities import cythonized
+
def poly_LC(f, K):
"""Returns leading coefficient of `f`. """
if not f:
@@ -24,13 +26,15 @@ def poly_TC(f, K):
dup_LC = dmp_LC = poly_LC
dup_TC = dmp_TC = poly_TC
+@cythonized("u")
def dmp_ground_LC(f, u, K):
"""Returns ground leading coefficient. """
if not u:
return dup_LC(f, K)
else:
return dmp_ground_LC(dmp_LC(f, K), u-1, K)
+@cythonized("u")
def dmp_ground_TC(f, u, K):
"""Returns ground trailing coefficient. """
if not u:
@@ -42,42 +46,53 @@ def dup_degree(f):
"""Returns leading degree of `f` in `K[x]`. """
return len(f) - 1
+@cythonized("u")
def dmp_degree(f, u):
"""Returns leading degree of `f` in `x_0` in `K[X]`. """
if dmp_zero_p(f, u):
return -1
else:
return len(f) - 1
+@cythonized("v,i,j")
+def _rec_degree_in(g, v, i, j):
+ """XXX"""
+ if i == j:
+ return dmp_degree(g, v)
+
+ v, i = v-1, i+1
+
+ return max([ _rec_degree_in(c, v, i, j) for c in g ])
+
+@cythonized("j,u")
def dmp_degree_in(f, j, u):
"""Returns leading degree of `f` in `x_j` in `K[X]`. """
if not j:
return dmp_degree(f, u)
if j < 0 or j > u:
raise IndexError("-%s <= j < %s expected, got %s" % (u, u, j))
- def rec_degree_in(g, v, i):
- if i == j:
- return dmp_degree(g, v)
- else:
- return max([ rec_degree_in(c, v-1, i+1) for c in g ])
+ return _rec_degree_in(f, u, 0, j)
- return rec_degree_in(f, u, 0)
+@cythonized("v,i")
+def _rec_degree_list(g, v, i, degs):
+ """XXX"""
+ degs[i] = max(degs[i], dmp_degree(g, v))
+ if v > 0:
+ v, i = v-1, i+1
+
+ for c in g:
+ _rec_degree_list(c, v, i, degs)
+
+@cythonized("u")
def dmp_degree_list(f, u):
"""Returns a list of degrees of `f` in `K[X]`. """
degs = [-1]*(u+1)
-
- def rec_degree_list(g, v, i):
- degs[i] = max(degs[i], dmp_degree(g, v))
-
- if v > 0:
- for c in g:
- rec_degree_list(c, v-1, i+1)
-
- rec_degree_list(f, u, 0)
+ _rec_degree_list(f, u, 0, degs)
return tuple(degs)
+@cythonized("i")
def dup_strip(f):
"""Remove leading zeros from `f` in `K[x]`. """
if not f or f[0]:
@@ -93,6 +108,7 @@ def dup_strip(f):
return f[i:]
+@cythonized("u,v,i")
def dmp_strip(f, u):
"""Remove leading zeros from `f` in `K[X]`. """
if not u:
@@ -114,57 +130,73 @@ def dmp_strip(f, u):
else:
return f[i:]
-def dmp_validate(f, K=None):
- """Returns number of levels in `f` and recursively strips it. """
- levels = set([])
+@cythonized("i,j")
+def _rec_validate(f, g, i, K):
+ """XXX"""
+ if type(g) is not list:
+ if K is not None and not K.of_type(g):
+ raise TypeError("%s in %s in not of type %s" % (g, f, K.dtype))
- def rec_validate(g, i):
- if type(g) is not list:
- if K is not None and not K.of_type(g):
- raise TypeError("%s in %s in not of type %s" % (g, f, K.dtype))
+ return set([i-1])
+ elif not g:
+ return set([i])
+ else:
+ j, levels = i+1, set([])
- levels.add(i-1)
- elif not g:
- levels.add(i)
- else:
- for c in g:
- rec_validate(c, i+1)
+ for c in g:
+ levels |= _rec_validate(f, c, i+1, K)
- def rec_strip(g, v):
- if not v:
- return dup_strip(g)
- else:
- return dmp_strip([ rec_strip(c, v-1) for c in g ], v)
+ return levels
+
+@cythonized("v,w")
+def _rec_strip(g, v):
+ """XXX"""
+ if not v:
+ return dup_strip(g)
+
+ w = v-1
+
+ return dmp_strip([ _rec_strip(c, w) for c in g ], v)
+
+@cythonized("u")
+def dmp_validate(f, K=None):
+ """Returns number of levels in `f` and recursively strips it. """
+ levels = _rec_validate(f, f, 0, K)
- rec_validate(f, 0)
u = levels.pop()
if not levels:
- return rec_strip(f, u), u
+ return _rec_strip(f, u), u
else:
raise ValueError("invalid data structure for a multivariate polynomial")
def dup_copy(f):
"""Create a new copy of a polynomial `f` in `K[x]`. """
return list(f)
+@cythonized("u,v")
def dmp_copy(f, u):
"""Create a new copy of a polynomial `f` in `K[X]`. """
if not u:
return list(f)
- else:
- return [ dmp_copy(c, u-1) for c in f ]
+
+ v = u-1
+
+ return [ dmp_copy(c, v) for c in f ]
def dup_normal(f, K):
"""Normalize univariate polynomial in the given domain. """
return dup_strip([ K(c) for c in f ])
+@cythonized("u,v")
def dmp_normal(f, u, K):
"""Normalize multivariate polynomial in the given domain. """
if not u:
return dup_normal(f, K)
- else:
- return dmp_strip([ dmp_normal(c, u-1, K) for c in f ], u)
+
+ v = u-1
+
+ return dmp_strip([ dmp_normal(c, v, K) for c in f ], u)
def dup_convert(f, K0, K1):
"""Convert ground domain of `f` from `K0` to `K1`. """
@@ -173,15 +205,19 @@ def dup_convert(f, K0, K1):
else:
return dup_strip([ K1.convert(c, K0) for c in f ])
+@cythonized("u,v")
def dmp_convert(f, u, K0, K1):
"""Convert ground domain of `f` from `K0` to `K1`. """
if K0 == K1:
return f
- elif not u:
+ if not u:
return dup_convert(f, K0, K1)
- else:
- return dmp_strip([ dmp_convert(c, u-1, K0, K1) for c in f ], u)
+ v = u-1
+
+ return dmp_strip([ dmp_convert(c, v, K0, K1) for c in f ], u)
+
+@cythonized("n")
def dup_nth(f, n, K):
"""Returns n-th coefficient of `f` in `K[x]`. """
if n < 0:
@@ -191,6 +227,7 @@ def dup_nth(f, n, K):
else:
return f[dup_degree(f)-n]
+@cythonized("n,u")
def dmp_nth(f, n, u, K):
"""Returns n-th coefficient of `f` in `K[x]`. """
if n < 0:
@@ -200,6 +237,7 @@ def dmp_nth(f, n, u, K):
else:
return f[dmp_degree(f, u)-n]
+@cythonized("n,u,v")
def dmp_ground_nth(f, N, u, K):
"""Returns ground n-th coefficient of `f` in `K[x]`. """
v = u
@@ -214,6 +252,7 @@ def dmp_ground_nth(f, N, u, K):
return f
+@cythonized("u")
def dmp_zero_p(f, u):
"""Returns True if `f` is zero in `K[X]`. """
if not u:
@@ -224,21 +263,25 @@ def dmp_zero_p(f, u):
else:
return False
+@cythonized("u")
def dmp_zero(u):
"""Returns a multivariate zero. """
if not u:
return []
else:
return [dmp_zero(u-1)]
+@cythonized("u")
def dmp_one_p(f, u, K):
"""Returns True if `f` is one in `K[X]`. """
return dmp_ground_p(f, K.one, u)
+@cythonized("u")
def dmp_one(u, K):
"""Returns a multivariate one over `K`. """
return dmp_ground(K.one, u)
+@cythonized("u")
def dmp_ground_p(f, c, u):
"""Returns True if `f` is constant in `K[X]`. """
if c is not None:
@@ -255,6 +298,7 @@ def dmp_ground_p(f, c, u):
else:
return False
+@cythonized("u")
def dmp_ground(c, u):
"""Returns a multivariate constant. """
if not c:
@@ -267,6 +311,7 @@ def dmp_ground(c, u):
else:
return [dmp_ground(c, u-1)]
+@cythonized("n,u")
def dmp_zeros(n, u, K):
"""Returns a list of multivariate zeros. """
if not n:
@@ -277,6 +322,7 @@ def dmp_zeros(n, u, K):
else:
return [ dmp_zero(u) for i in xrange(n) ]
+@cythonized("n,u")
def dmp_grounds(c, n, u):
"""Returns a list of multivariate constants. """
if not n:
@@ -287,26 +333,30 @@ def dmp_grounds(c, n, u):
else:
return [ dmp_ground(c, u) for i in xrange(n) ]
+@cythonized("u")
def dmp_negative_p(f, u, K):
"""Returns `True` if `LC(f)` is negative. """
return K.is_negative(dmp_ground_LC(f, u, K))
+@cythonized("u")
def dmp_positive_p(f, u, K):
"""Returns `True` if `LC(f)` is positive. """
return K.is_positive(dmp_ground_LC(f, u, K))
+@cythonized("n,k")
def dup_from_dict(f, K):
"""Create `K[x]` polynomial from a dict. """
if not f:
return []
- n, h = max(k for (k,) in f.iterkeys()), []
+ n, h = max([ k for (k,) in f.iterkeys() ]), []
for k in xrange(n, -1, -1):
h.append(f.get((k,), K.zero))
return dup_strip(h)
+@cythonized("n,k")
def dup_from_raw_dict(f, K):
"""Create `K[x]` polynomial from a raw dict. """
if not f:
@@ -319,6 +369,7 @@ def dup_from_raw_dict(f, K):
return dup_strip(h)
+@cythonized("u,v,n,k")
def dmp_from_dict(f, u, K):
"""Create `K[X]` polynomial from a dict. """
if not u:
@@ -348,41 +399,45 @@ def dmp_from_dict(f, u, K):
return dmp_strip(h, u)
+@cythonized("n,k")
def dup_to_dict(f):
"""Convert `K[x]` polynomial to a dict. """
n, result = dup_degree(f), {}
- for i in xrange(0, n+1):
- if f[n-i]:
- result[(i,)] = f[n-i]
+ for k in xrange(0, n+1):
+ if f[n-k]:
+ result[(k,)] = f[n-k]
return result
+@cythonized("n,k")
def dup_to_raw_dict(f):
"""Convert `K[x]` polynomial to a raw dict. """
n, result = dup_degree(f), {}
- for i in xrange(0, n+1):
- if f[n-i]:
- result[i] = f[n-i]
+ for k in xrange(0, n+1):
+ if f[n-k]:
+ result[k] = f[n-k]
return result
+@cythonized("u,v,n,k")
def dmp_to_dict(f, u):
"""Convert `K[X]` polynomial to a dict. """
if not u:
return dup_to_dict(f)
n, v, result = dmp_degree(f, u), u-1, {}
- for i in xrange(0, n+1):
- h = dmp_to_dict(f[n-i], v)
+ for k in xrange(0, n+1):
+ h = dmp_to_dict(f[n-k], v)
for exp, coeff in h.iteritems():
- result[(i,)+exp] = coeff
+ result[(k,)+exp] = coeff
return result
+@cythonized("u,i,j")
def dmp_swap(f, i, j, u, K):
"""Transform `K[..x_i..x_j..]` to `K[..x_j..x_i..]`. """
if i < 0 or j < 0 or i > u or j > u:
@@ -399,6 +454,7 @@ def dmp_swap(f, i, j, u, K):
return dmp_from_dict(H, u, K)
+@cythonized("u")
def dmp_permute(f, P, u, K):
"""Returns a polynomial in `K[x_{P(1)},..,x_{P(n)}]`. """
F, H = dmp_to_dict(f, u), {}
@@ -413,6 +469,7 @@ def dmp_permute(f, P, u, K):
return dmp_from_dict(H, u, K)
+@cythonized("l")
def dmp_nest(f, l, K):
"""Returns multivariate value nested `l`-levels. """
if type(f) is not list:
@@ -423,6 +480,7 @@ def dmp_nest(f, l, K):
else:
return [dmp_nest(f, l-1, K)]
+@cythonized("l,k,u,v")
def dmp_raise(f, l, u, K):
"""Returns multivariate polynomial raised `l`-levels. """
if not l:
@@ -431,11 +489,16 @@ def dmp_raise(f, l, u, K):
if not u:
if not f:
return dmp_zero(l)
- else:
- return [ dmp_ground(c, l-1) for c in f ]
- else:
- return [ dmp_raise(c, l, u-1, K) for c in f ]
+ k = l-1
+
+ return [ dmp_ground(c, l-1) for c in f ]
+
+ v = u-1
+
+ return [ dmp_raise(c, l, u-1, K) for c in f ]
+
+@cythonized("g,i")
def dup_deflate(f, K):
"""Map `x**m` to `y` in a polynomial in `K[x]`. """
if dup_degree(f) <= 0:
@@ -454,35 +517,37 @@ def dup_deflate(f, K):
return g, f[::g]
+@cythonized("u,i,m,a,b")
def dmp_deflate(f, u, K):
"""Map `x_i**m_i` to `y_i` in a polynomial in `K[X]`. """
if dmp_zero_p(f, u):
return (1,)*(u+1), f
- F, H = dmp_to_dict(f, u), {}
-
- def ilgcd(M):
- g = 0
+ F = dmp_to_dict(f, u)
+ B = [0]*(u+1)
- for m in M:
- g = igcd(g, m)
+ for M in F.iterkeys():
+ for i, m in enumerate(M):
+ B[i] = igcd(B[i], m)
- if g == 1:
- break
+ for i, b in enumerate(B):
+ if not b:
+ B[i] = 1
- return g or 1
+ B = tuple(B)
- M = tuple(map(lambda *row: ilgcd(row), *F.keys()))
+ if all([ b == 1 for b in B ]):
+ return B, f
- if all([ b == 1 for b in M ]):
- return M, f
+ H = {}
- for m, coeff in F.iteritems():
- N = [ a // b for a, b in zip(m, M) ]
+ for A, coeff in F.iteritems():
+ N = [ a // b for a, b in zip(A, B) ]
H[tuple(N)] = coeff
- return M, dmp_from_dict(H, u, K)
+ return B, dmp_from_dict(H, u, K)
+@cythonized("G,g,i")
def dup_multi_deflate(polys, K):
"""Map `x**m` to `y` in a set of polynomials in `K[x]`. """
G = 0
@@ -506,52 +571,48 @@ def dup_multi_deflate(polys, K):
return G, tuple([ p[::G] for p in polys ])
+@cythonized("u,G,g,m,a,b")
def dmp_multi_deflate(polys, u, K):
"""Map `x_i**m_i` to `y_i` in a set of polynomials in `K[X]`. """
- def ilgcd(M):
- g = 0
-
- for m in M:
- g = igcd(g, m)
-
- if g == 1:
- break
-
- return g or 1
-
if not u:
M, H = dup_multi_deflate(polys, K)
return (M,), H
- F, M, H = [], [], []
+ F, B = [], [0]*(u+1)
for p in polys:
f = dmp_to_dict(p, u)
- if dmp_zero_p(p, u):
- m = (0,)*(u+1)
- else:
- m = map(lambda *row: ilgcd(row), *f.keys())
+ if not dmp_zero_p(p, u):
+ for M in f.iterkeys():
+ for i, m in enumerate(M):
+ B[i] = igcd(B[i], m)
F.append(f)
- M.append(m)
- M = tuple(map(lambda *row: ilgcd(row), *M))
+ for i, b in enumerate(B):
+ if not b:
+ B[i] = 1
+
+ B = tuple(B)
- if all([ b == 1 for b in M ]):
- return M, polys
+ if all([ b == 1 for b in B ]):
+ return B, polys
+
+ H = []
for f in F:
h = {}
- for m, coeff in f.iteritems():
- N = [ a // b for a, b in zip(m, M) ]
+ for A, coeff in f.iteritems():
+ N = [ a // b for a, b in zip(A, B) ]
h[tuple(N)] = coeff
H.append(dmp_from_dict(h, u, K))
- return M, tuple(H)
+ return B, tuple(H)
+@cythonized("m")
def dup_inflate(f, m, K):
"""Map `y` to `x**m` in a polynomial in `K[x]`. """
if m <= 0:
@@ -567,34 +628,40 @@ def dup_inflate(f, m, K):
return result
-def dmp_inflate(f, M, u, K):
- """Map `y_i` to `x_i**k_i` in a polynomial in `K[X]`. """
- if not u:
- return dup_inflate(f, M[0], K)
+@cythonized("u,v,i,j")
+def _rec_inflate(g, M, v, i, K):
+ """XXX"""
+ if not v:
+ return dup_inflate(g, M[i], K)
+ if M[i] <= 0:
+ raise IndexError("all `M[i]` must be positive, got %s" % M[i])
- def rec_inflate(g, v, i):
- if not v:
- return dup_inflate(g, M[i], K)
+ w, j = v-1, i+1
- if M[i] <= 0:
- raise IndexError("all `M[i]` must be positive, got %s" % M[i])
+ g = [ _rec_inflate(c, M, w, j, K) for c in g ]
- g = [ rec_inflate(c, v-1, i+1) for c in g ]
- result, w = [g[0]], v-1
+ result = [g[0]]
- for coeff in g[1:]:
- for _ in xrange(1, M[i]):
- result.append(dmp_zero(w))
+ for coeff in g[1:]:
+ for _ in xrange(1, M[i]):
+ result.append(dmp_zero(w))
- result.append(coeff)
+ result.append(coeff)
- return result
+ return result
+
+@cythonized("u,m")
+def dmp_inflate(f, M, u, K):
+ """Map `y_i` to `x_i**k_i` in a polynomial in `K[X]`. """
+ if not u:
+ return dup_inflate(f, M[0], K)
if all([ m == 1 for m in M ]):
return f
else:
- return rec_inflate(f, u, 0)
+ return _rec_inflate(f, M, u, 0, K)
+@cythonized("u,j")
def dmp_exclude(f, u, K):
"""Exclude useless levels from `f`. """
if not u or dmp_ground_p(f, None, u):
@@ -626,6 +693,7 @@ def dmp_exclude(f, u, K):
return J, dmp_from_dict(f, u, K), u
+@cythonized("u,j")
def dmp_include(f, J, u, K):
"""Include useless levels in `f`. """
if not J:
@@ -645,6 +713,7 @@ def dmp_include(f, J, u, K):
return dmp_from_dict(f, u, K)
+@cythonized("u,v,w")
def dmp_inject(f, u, K, **args):
"""Convert `f` from `K[X][Y]` to `K[X,Y]`. """
front = args.get('front', False)
@@ -665,6 +734,7 @@ def dmp_inject(f, u, K, **args):
return dmp_from_dict(h, w, K.dom), w
+@cythonized("u,v")
def dmp_eject(f, u, K, **args):
"""Convert `f` from `K[X,Y]` to `K[X][Y]`. """
front = args.get('front', False)
@@ -688,68 +758,74 @@ def dmp_eject(f, u, K, **args):
return dmp_from_dict(h, v-1, K)
+@cythonized("g,G")
def dup_terms_gcd(f, K):
"""Remove GCD of terms from `f` in `K[x]`. """
if dup_TC(f, K) or not f:
return 0, f
F = dup_to_raw_dict(f)
+ G = min(F.keys())
- gcd = min(F.keys())
-
- if not gcd:
- return gcd, f
+ if not G:
+ return G, f
f = {}
- for k, coeff in F.iteritems():
- f[k - gcd] = coeff
+ for g, coeff in F.iteritems():
+ f[g-G] = coeff
- return gcd, dup_from_raw_dict(f, K)
+ return G, dup_from_raw_dict(f, K)
+@cythonized("u,g")
def dmp_terms_gcd(f, u, K):
"""Remove GCD of terms from `f` in `K[X]`. """
if dmp_ground_TC(f, u, K) or dmp_zero_p(f, u):
return (0,)*(u+1), f
F = dmp_to_dict(f, u)
+ G = monomial_min(*F.keys())
- gcd = monomial_min(*F.keys())
-
- if all(not n for n in gcd):
- return gcd, f
+ if all([ g == 0 for g in G ]):
+ return G, f
f = {}
for monom, coeff in F.iteritems():
- f[monomial_div(monom, gcd)] = coeff
+ f[monomial_div(monom, G)] = coeff
- return gcd, dmp_from_dict(f, u, K)
+ return G, dmp_from_dict(f, u, K)
-def dmp_list_terms(f, u, K):
- """List all non-zero terms from `f` in lex order. """
- terms = []
+@cythonized("v,w,d,i")
+def _rec_list_terms(g, v, monom):
+ """XXX"""
+ d, terms = dmp_degree(g, v), []
- def rec_iter_terms(g, v, monom):
- d = dmp_degree(g, v)
+ if not v:
+ for i, c in enumerate(g):
+ if not c:
+ continue
- if not v:
- for i, c in enumerate(g):
- if not c:
- continue
+ terms.append((monom + (d-i,), c))
+ else:
+ w = v-1
- terms.append((monom + (d-i,), c))
- else:
- for i, c in enumerate(g):
- rec_iter_terms(c, v-1, monom + (d-i,))
+ for i, c in enumerate(g):
+ terms.extend(_rec_list_terms(c, v-1, monom + (d-i,)))
+
+ return terms
- rec_iter_terms(f, u, ())
+@cythonized("u")
+def dmp_list_terms(f, u, K):
+ """List all non-zero terms from `f` in lex order. """
+ terms = _rec_list_terms(f, u, ())
if not terms:
return [((0,)*(u+1), K.zero)]
else:
return terms
+@cythonized("n,m")
def dup_apply_pairs(f, g, h, args, K):
"""Apply `h` to pairs of coefficients of `f` and `g`. """
n, m = len(f), len(g)
@@ -767,6 +843,7 @@ def dup_apply_pairs(f, g, h, args, K):
return dup_strip(result)
+@cythonized("u,v,n,m")
def dmp_apply_pairs(f, g, h, args, u, K):
"""Apply `h` to pairs of coefficients of `f` and `g`. """
if not u:
View
452 sympy/polys/densetools.py
@@ -11,7 +11,7 @@
dmp_multi_deflate, dmp_inflate,
dup_to_raw_dict, dup_from_raw_dict,
dmp_raise, dmp_apply_pairs,
- dmp_zeros,
+ dmp_inject, dmp_zeros
)
from sympy.polys.densearith import (
@@ -31,21 +31,22 @@
dup_mul_ground, dmp_mul_ground,
dup_quo_ground, dmp_quo_ground,
dup_exquo_ground, dmp_exquo_ground,
- dup_max_norm, dmp_max_norm,
+ dup_max_norm, dmp_max_norm
)
from sympy.polys.galoistools import (
- gf_int, gf_crt,
+ gf_int, gf_crt
)
from sympy.polys.polyerrors import (
HeuristicGCDFailed,
HomomorphismFailed,
NotInvertible,
- DomainError,
+ DomainError
)
from sympy.ntheory import nextprime
+from sympy.utilities import cythonized
def dup_ground_to_ring(f, K0, K1):
"""Clear denominators, i.e. transform `K_0` to `K_1`, but don't normalize. """
@@ -59,30 +60,36 @@ def dup_ground_to_ring(f, K0, K1):
else:
return common, dup_mul_ground(f, common, K0)
+@cythonized("v,w")
+def _rec_ground_to_ring(g, v, K0, K1):
+ """XXX"""
+ common = K1.one
+
+ if not v:
+ for c in g:
+ common = K1.lcm(common, K0.denom(c))
+ else:
+ w = v-1
+
+ for c in g:
+ common = K1.lcm(common, _rec_ground_to_ring(c, w, K0, K1))
+
+ return common
+
+@cythonized("u")
def dmp_ground_to_ring(f, u, K0, K1):
"""Clear denominators, i.e. transform `K_0` to `K_1`, but don't normalize. """
if not u:
return dup_ground_to_ring(f, K0, K1)
- def rec_ground_to_ring(g, v):
- common = K1.one
-
- if not v:
- for c in g:
- common = K1.lcm(common, K0.denom(c))
- else:
- for c in g:
- common = K1.lcm(common, rec_ground_to_ring(c, v-1))
-
- return common
-
- common = rec_ground_to_ring(f, u)
+ common = _rec_ground_to_ring(f, u, K0, K1)
if not K1.is_one(common):
f = dmp_mul_ground(f, common, u, K0)
return common, f
+@cythonized("m,n,i,j")
def dup_integrate(f, m, K):
"""Computes indefinite integral of `f` in `K[x]`. """
if m <= 0 or not f:
@@ -96,10 +103,11 @@ def dup_integrate(f, m, K):
for j in xrange(1, m):
n *= i+j+1
- g.insert(0, K.quo(c, n))
+ g.insert(0, K.quo(c, K(n)))
return g
+@cythonized("m,u,v,n,i,j")
def dmp_integrate(f, m, u, K):
"""Computes indefinite integral of `f` in `x_0` in `K[X]`. """
if not u:
@@ -108,31 +116,37 @@ def dmp_integrate(f, m, u, K):
if m <= 0 or dmp_zero_p(f, u):
return f
- g = dmp_zeros(m, u-1, K)
+ g, v = dmp_zeros(m, u-1, K), u-1
for i, c in enumerate(reversed(f)):
n = i+1
for j in xrange(1, m):
n *= i+j+1
- g.insert(0, dmp_quo_ground(c, K(n), u-1, K))
+ g.insert(0, dmp_quo_ground(c, K(n), v, K))
return g
+@cythonized("m,v,w,i,j")
+def _rec_integrate_in(g, m, v, i, j, K):
+ """XXX"""
+ if i == j:
+ return dmp_integrate(g, m, v, K)
+
+ w, i = v-1, i+1
+
+ return dmp_strip([ _rec_integrate_in(c, m, w, i, j, K) for c in g ], v)
+
+@cythonized("m,j,u")
def dmp_integrate_in(f, m, j, u, K):
"""Computes indefinite integral of `f` in `x_j` in `K[X]`. """
if j < 0 or j > u:
raise IndexError("-%s <= j < %s expected, got %s" % (u, u, j))
- def rec_integrate_in(g, v, i):
- if i == j:
- return dmp_integrate(g, m, v, K)
- else:
- return dmp_strip([ rec_integrate_in(c, v-1, i+1) for c in g ], v)
-
- return rec_integrate_in(f, u, 0)
+ return _rec_integrate_in(f, m, u, 0, j, K)
+@cythonized("m,n,i")
def dup_diff(f, m, K):
"""m-th order derivative of a polynomial in `K[x]`. """
if m <= 0:
@@ -146,14 +160,15 @@ def dup_diff(f, m, K):
deriv, c = [], K.one
for i in xrange(0, m):
- c, n = c*n, n-1
+ c, n = c*K(n), n-1
for coeff in f[:-m]:
deriv.append(coeff*c)
- c, n = n*K.exquo(c, n+m), n-1
+ c, n = K(n)*K.exquo(c, K(n+m)), n-1
return deriv
+@cythonized("u,v,m,n,i")
def dmp_diff(f, m, u, K):
"""m-th order derivative in `x_0` of a polynomial in `K[X]`. """
if not u:
@@ -169,27 +184,32 @@ def dmp_diff(f, m, u, K):
deriv, c, v = [], K.one, u-1
for i in xrange(0, m):
- c, n = c*n, n-1
+ c, n = c*K(n), n-1
for coeff in f[:-m]:
h = dmp_mul_ground(coeff, c, v, K)
- c, n = n*K.exquo(c, n+m), n-1
+ c, n = K(n)*K.exquo(c, K(n+m)), n-1
deriv.append(h)
return deriv
+@cythonized("m,v,w,i,j")
+def _rec_diff_in(g, m, v, i, j, K):
+ """XXX"""
+ if i == j:
+ return dmp_diff(g, m, v, K)
+
+ w, i = v-1, i+1
+
+ return dmp_strip([ _rec_diff_in(c, m, w, i, j, K) for c in g ], v)
+
+@cythonized("m,j,u")
def dmp_diff_in(f, m, j, u, K):
"""m-th order derivative in `x_j` of a polynomial in `K[X]`. """
if j < 0 or j > u:
raise IndexError("-%s <= j < %s expected, got %s" % (u, u, j))
- def rec_diff_in(g, v, i):
- if i == j:
- return dmp_diff(g, m, v, K)
- else:
- return dmp_strip([ rec_diff_in(c, v-1, i+1) for c in g ], v)
-
- return rec_diff_in(f, u, 0)
+ return _rec_diff_in(f, m, u, 0, j, K)
def dup_eval(f, a, K):
"""Evaluate a polynomial at `x = a` in `K[x]` using Horner scheme. """
@@ -204,6 +224,7 @@ def dup_eval(f, a, K):
return result
+@cythonized("u,v")
def dmp_eval(f, a, u, K):
"""Evaluate a polynomial at `x_0 = a` in `K[X]` using Horner scheme. """
if not u:
@@ -220,62 +241,72 @@ def dmp_eval(f, a, u, K):
return result
+@cythonized("v,i,j")
+def _rec_eval_in(g, a, v, i, j, K):
+ """XXX"""
+ if i == j:
+ return dmp_eval(g, a, v, K)
+
+ v, i = v-1, i+1
+
+ return dmp_strip([ _rec_eval_in(c, a, v, i, j, K) for c in g ], v)
+
+@cythonized("u")
def dmp_eval_in(f, a, j, u, K):
"""Evaluate a polynomial at `x_j = a` in `K[X]` using Horner scheme. """
if j < 0 or j > u:
raise IndexError("-%s <= j < %s expected, got %s" % (u, u, j))
- def rec_eval_in(g, v, i):
- if i == j:
- return dmp_eval(g, a, v, K)
- else:
- return dmp_strip([ rec_eval_in(c, v-1, i+1) for c in g ], v-1)
+ return _rec_eval_in(f, a, u, 0, j, K)
- return rec_eval_in(f, u, 0)
+@cythonized("i,u")
+def _rec_eval_tail(g, i, A, u, K):
+ """XXX"""
+ if i == u:
+ return dup_eval(g, A[-1], K)
+ else:
+ h = [ _rec_eval_tail(c, i+1, A, u, K) for c in g ]
-def dmp_eval_tail(f, A, u, K):
- """Evaluate a polynomial at `x_j = a_j, ...` in `K[X]`. """
- def rec_eval_tail(g, i):
- if i == u:
- return dup_eval(g, A[-1], K)
+ if i < u - len(A) + 1:
+ return h
else:
- h = [ rec_eval_tail(c, i+1) for c in g ]
-
- if i < u - len(A) + 1:
- return h
- else:
- return dup_eval(h, A[-u+i-1], K)
+ return dup_eval(h, A[-u+i-1], K)
+@cythonized("u")
+def dmp_eval_tail(f, A, u, K):
+ """Evaluate a polynomial at `x_j = a_j, ...` in `K[X]`. """
if not A:
return f
if dmp_zero_p(f, u):
return dmp_zero(u - len(A))
- e = rec_eval_tail(f, 0)
+ e = _rec_eval_tail(f, 0, A, u, K)
if u == len(A)-1:
return e
else:
return dmp_strip(e, u - len(A))
+@cythonized("m,v,i,j")
+def _rec_diff_eval(g, m, a, v, i, j, K):
+ """XXX"""
+ if i == j:
+ return dmp_eval(dmp_diff(g, m, v, K), a, v, K)
+
+ v, i = v-1, i+1
+
+ return dmp_strip([ _rec_diff_eval(c, m, a, v, i, j, K) for c in g ], v)
+
+@cythonized("m,j,u")
def dmp_diff_eval_in(f, m, a, j, u, K):
"""Differentiate and evaluate a polynomial in `x_j` at `a` in `K[X]`. """
- if j < 0:
- j += u + 1
-
if j > u:
raise IndexError("-%s <= j < %s expected, got %s" % (u, u, j))
if not j:
return dmp_eval(dmp_diff(f, m, u, K), a, u, K)
- def rec_diff_eval(g, v, i):
- if i == j:
- return dmp_eval(dmp_diff(g, m, v, K), a, v, K)
- else:
- return dmp_strip([ rec_diff_eval(c, v-1, i+1) for c in g ], v-1)
-
- return rec_diff_eval(f, u, 0)
+ return _rec_diff_eval(f, m, a, u, 0, j, K)
def dup_half_gcdex(f, g, K):
"""Half extended Euclidean algorithm in `F[x]`. """
@@ -312,6 +343,7 @@ def dup_invert(f, g, K):
else:
raise NotInvertible("zero divisor")
+@cythonized("n,m,d,k")
def dup_inner_subresultants(f, g, K):
"""Subresultant PRS algorithm in `K[x]`. """
n = dup_degree(f)
@@ -363,6 +395,7 @@ def dup_subresultants(f, g, K):
"""Computes subresultant PRS of two polynomials in `K[x]`. """
return dup_inner_subresultants(f, g, K)[0]
+@cythonized("s,i,du,dv,dw")
def dup_prs_resultant(f, g, K):
"""Resultant algorithm in `K[x]` using subresultant PRS. """
if not f or not g:
@@ -394,9 +427,9 @@ def dup_prs_resultant(f, g, K):
if s < 0:
p = -p
- j = dup_degree(R[-2])
+ i = dup_degree(R[-2])
- res = dup_LC(R[-1], K)**j
+ res = dup_LC(R[-1], K)**i
res = K.quo(res*p, q)
return res, R
@@ -405,6 +438,7 @@ def dup_resultant(f, g, K):
"""Computes resultant of two polynomials in `K[x]`. """
return dup_prs_resultant(f, g, K)[0]
+@cythonized("u,v,n,m,d,k")
def dmp_inner_subresultants(f, g, u, K):
"""Subresultant PRS algorithm in `K[X]`. """
if not u:
@@ -459,10 +493,12 @@ def dmp_inner_subresultants(f, g, u, K):
return R, B, D
+@cythonized("u")
def dmp_subresultants(f, g, u, K):
"""Computes subresultant PRS of two polynomials in `K[X]`. """
return dmp_inner_subresultants(f, g, u, K)[0]
+@cythonized("u,v,s,i,d,du,dv,dw")
def dmp_prs_resultant(f, g, u, K):
"""Resultant algorithm in `K[X]` using subresultant PRS. """
if not u:
@@ -502,13 +538,14 @@ def dmp_prs_resultant(f, g, u, K):
if s < 0:
p = dmp_neg(p, v, K)
- j = dmp_degree(R[-2], u)
+ i = dmp_degree(R[-2], u)
- res = dmp_pow(dmp_LC(R[-1], K), j, v, K)
+ res = dmp_pow(dmp_LC(R[-1], K), i, v, K)
res = dmp_quo(dmp_mul(res, p, v, K), q, v, K)
return res, R
+@cythonized("u,v,n,m,N,M,B")
def dmp_zz_modular_resultant(f, g, p, u, K):
"""Compute resultant of `f` and `g` modulo a prime `p`. """
if not u:
@@ -566,10 +603,13 @@ def dmp_zz_modular_resultant(f, g, p, u, K):
return r
+def _collins_crt(r, R, P, p, K):
+ """Wrapper of CRT for Collins's resultant algorithm. """
+ return gf_int(gf_crt([r, R], [P, p], K), P*p)
+
+@cythonized("u,v,n,m")
def dmp_zz_collins_resultant(f, g, u, K):
"""Collins's modular resultant algorithm in `Z[X]`. """
- def crt(r, R, P, p):
- return gf_int(gf_crt([r, R], [P, p], K), P*p)
n = dmp_degree(f, u)
m = dmp_degree(g, u)
@@ -605,12 +645,13 @@ def crt(r, R, P, p):
if K.is_one(P):
r = R
else:
- r = dmp_apply_pairs(r, R, crt, (P, p), v, K)
+ r = dmp_apply_pairs(r, R, _collins_crt, (P, p, K), v, K)
P *= p
return r
+@cythonized("u,n,m")
def dmp_qq_collins_resultant(f, g, u, K0):
"""Collins's modular resultant algorithm in `Q[X]`. """
n = dmp_degree(f, u)
@@ -636,6 +677,7 @@ def dmp_qq_collins_resultant(f, g, u, K0):
USE_COLLINS_RESULTANT = 0
+@cythonized("u")
def dmp_resultant(f, g, u, K):
"""Computes resultant of two polynomials in `K[X]`. """
if not u:
@@ -650,6 +692,7 @@ def dmp_resultant(f, g, u, K):
return dmp_prs_resultant(f, g, u, K)[0]
+@cythonized("d,s")
def dup_discriminant(f, K):
"""Computes discriminant of a polynomial in `K[x]`. """
d = dup_degree(f)
@@ -662,8 +705,9 @@ def dup_discriminant(f, K):
r = dup_resultant(f, dup_diff(f, 1, K), K)
- return K.quo(r, c*s)
+ return K.quo(r, c*K(s))
+@cythonized("u,v,d,s")
def dmp_discriminant(f, u, K):
"""Computes discriminant of a polynomial in `K[X]`. """
if not u:
@@ -678,7 +722,7 @@ def dmp_discriminant(f, u, K):
c = dmp_LC(f, K)
r = dmp_resultant(f, dmp_diff(f, 1, u, K), u, K)
- c = dmp_mul_ground(c, s, v, K)
+ c = dmp_mul_ground(c, K(s), v, K)
return dmp_quo(r, c, v, K)
@@ -712,6 +756,7 @@ def _dup_ff_trivial_gcd(f, g, K):
USE_DMP_SIMPLIFY_GCD = 1
+@cythonized("u")
def _dmp_rr_trivial_gcd(f, g, u, K):
"""Handle trivial cases in GCD algorithm over a ring. """
zero_f = dmp_zero_p(f, u)
@@ -734,6 +779,7 @@ def _dmp_rr_trivial_gcd(f, g, u, K):
else:
return None
+@cythonized("u")
def _dmp_ff_trivial_gcd(f, g, u, K):
"""Handle trivial cases in GCD algorithm over a field. """
zero_f = dmp_zero_p(f, u)
@@ -754,6 +800,7 @@ def _dmp_ff_trivial_gcd(f, g, u, K):
else:
return None
+@cythonized("u,v,df,dg")
def _dmp_simplify_gcd(f, g, u, K):
"""Try to eliminate `x_0` from GCD computation in `K[X]`. """
df = dmp_degree(f, u)
@@ -821,6 +868,7 @@ def dup_ff_prs_gcd(f, g, K):
return h, cff, cfg
+@cythonized("u")
def dmp_rr_prs_gcd(f, g, u, K):
"""Computes polynomial GCD using subresultants over a ring. """
if not u:
@@ -848,6 +896,7 @@ def dmp_rr_prs_gcd(f, g, u, K):
return h, cff, cfg
+@cythonized("u")
def dmp_ff_prs_gcd(f, g, u, K):
"""Computes polynomial GCD using subresultants over a field. """
if not u:
@@ -875,6 +924,22 @@ def dmp_ff_prs_gcd(f, g, u, K):
HEU_GCD_MAX = 6
+def _dup_zz_gcd_interpolate(h, x, K):
+ """Interpolate polynomial GCD from integer GCD. """
+ f = []
+
+ while h:
+ g = h % x
+
+ if g > x // 2:
+ g -= x
+
+ f.insert(0, g)
+ h = (h-g) // x
+
+ return f
+
+@cythonized("i,df,dg")
def dup_zz_heu_gcd(f, g, K):
"""Heuristic polynomial GCD in `Z[x]`.
@@ -907,24 +972,6 @@ def dup_zz_heu_gcd(f, g, K):
if result is not None:
return result
- def interpolate(h, x):
- f = []
-
- while h:
- g = h % x
-
- if g > x // 2:
- g -= x
-
- f.insert(0, g)
- h = (h-g) // x
-
- return f
-
- def finalize(h, cff, cfg, gcd):
- h = dup_mul_ground(h, gcd, K)
- return h, cff, cfg
-
df = dup_degree(f)
dg = dup_degree(g)
@@ -952,7 +999,7 @@ def finalize(h, cff, cfg, gcd):
cff = ff // h
cfg = gg // h
- h = interpolate(h, x)
+ h = _dup_zz_gcd_interpolate(h, x, K)
h = dup_primitive(h, K)[1]
cff_, r = dup_div(f, h, K)
@@ -961,32 +1008,53 @@ def finalize(h, cff, cfg, gcd):
cfg_, r = dup_div(g, h, K)
if not r:
- return finalize(h, cff_, cfg_, gcd)
+ h = dup_mul_ground(h, gcd, K)
+ return h, cff_, cfg_
- cff = interpolate(cff, x)
+ cff = _dup_zz_gcd_interpolate(cff, x, K)
h, r = dup_div(f, cff, K)
if not r:
cfg_, r = dup_div(g, h, K)
if not r:
- return finalize(h, cff, cfg_, gcd)
+ h = dup_mul_ground(h, gcd, K)
+ return h, cff, cfg_
- cfg = interpolate(cfg, x)
+ cfg = _dup_zz_gcd_interpolate(cfg, x, K)
h, r = dup_div(g, cfg, K)
if not r:
cff_, r = dup_div(f, h, K)
if not r:
- return finalize(h, cff_, cfg, gcd)
+ h = dup_mul_ground(h, gcd, K)
+ return h, cff, cfg
x = 73794*x * K.sqrt(K.sqrt(x)) // 27011
raise HeuristicGCDFailed('no luck')
+@cythonized("v")
+def _dmp_zz_gcd_interpolate(h, x, v, K):
+ """Interpolate polynomial GCD from integer GCD. """
+ f = []
+
+ while not dmp_zero_p(h, v):
+ g = dmp_ground_trunc(h, x, v, K)
+ f.insert(0, g)
+
+ h = dmp_sub(h, g, v, K)
+ h = dmp_exquo_ground(h, x, v, K)
+
+ if K.is_negative(dmp_ground_LC(f, v+1, K)):
+ return dmp_neg(f, v+1, K)
+ else:
+ return f
+
+@cythonized("u,v,i,dg,df")
def dmp_zz_heu_gcd(f, g, u, K):
"""Heuristic polynomial GCD in `Z[X]`.
@@ -1024,25 +1092,6 @@ def dmp_zz_heu_gcd(f, g, u, K):
if result is not None:
return result
- def interpolate(h, x, v):
- f = []
-
- while not dmp_zero_p(h, v):
- g = dmp_ground_trunc(h, x, v, K)
- f.insert(0, g)
-
- h = dmp_sub(h, g, v, K)
- h = dmp_exquo_ground(h, x, v, K)
-
- if K.is_negative(dmp_ground_LC(f, v+1, K)):
- return dmp_neg(f, v+1, K)
- else:
- return f
-
- def finalize(h, cff, cfg, gcd):
- h = dmp_mul_ground(h, gcd, u, K)
- return h, cff, cfg
-
df = dmp_degree(f, u)
dg = dmp_degree(g, u)
@@ -1066,7 +1115,7 @@ def finalize(h, cff, cfg, gcd):
if not (dmp_zero_p(ff, v) or dmp_zero_p(gg, v)):
h, cff, cfg = dmp_zz_heu_gcd(ff, gg, v, K)
- h = interpolate(h, x, v)
+ h = _dmp_zz_gcd_interpolate(h, x, v, K)
h = dmp_ground_primitive(h, u, K)[1]
cff_, r = dmp_div(f, h, u, K)
@@ -1075,27 +1124,30 @@ def finalize(h, cff, cfg, gcd):
cfg_, r = dmp_div(g, h, u, K)
if dmp_zero_p(r, u):
- return finalize(h, cff_, cfg_, gcd)
+ h = dmp_mul_ground(h, gcd, u, K)
+ return h, cff_, cfg_
- cff = interpolate(cff, x, v)
+ cff = _dmp_zz_gcd_interpolate(cff, x, v, K)
h, r = dmp_div(f, cff, u, K)
if dmp_zero_p(r, u):
cfg_, r = dmp_div(g, h, u, K)
if dmp_zero_p(r, u):
- return finalize(h, cff, cfg_, gcd)
+ h = dmp_mul_ground(h, gcd, u, K)
+ return h, cff, cfg_
- cfg = interpolate(cfg, x, v)
+ cfg = _dmp_zz_gcd_interpolate(cfg, x, v, K)
h, r = dmp_div(g, cfg, u, K)
if dmp_zero_p(r, u):
cff_, r = dmp_div(f, h, u, K)
if dmp_zero_p(r, u):
- return finalize(h, cff_, cfg, gcd)
+ h = dmp_mul_ground(h, gcd, u, K)
+ return h, cff_, cfg
x = 73794*x * K.sqrt(K.sqrt(x)) // 27011
@@ -1131,6 +1183,7 @@ def dup_qq_heu_gcd(f, g, K0):
return h, cff, cfg
+@cythonized("u")
def dmp_qq_heu_gcd(f, g, u, K0):
"""Heuristic polynomial GCD in `Q[X]`. """
result = _dmp_ff_trivial_gcd(f, g, u, K0)
@@ -1185,6 +1238,7 @@ def dup_inner_gcd(f, g, K):
return dup_rr_prs_gcd(f, g, K)
+@cythonized("u")
def _dmp_inner_gcd(f, g, u, K):
"""Helper function for `dmp_inner_gcd()`. """
if K.has_Field:
@@ -1206,6 +1260,7 @@ def _dmp_inner_gcd(f, g, u, K):
return dmp_rr_prs_gcd(f, g, u, K)
+@cythonized("u")
def dmp_inner_gcd(f, g, u, K):
"""Computes polynomial GCD and cofactors of `f` and `g` in `K[X]`. """
if not u:
@@ -1222,6 +1277,7 @@ def dup_gcd(f, g, K):
"""Computes polynomial GCD of `f` and `g` in `K[x]`. """
return dup_inner_gcd(f, g, K)[0]
+@cythonized("u")
def dmp_gcd(f, g, u, K):
"""Computes polynomial GCD of `f` and `g` in `K[X]`. """
return dmp_inner_gcd(f, g, u, K)[0]
@@ -1243,7 +1299,7 @@ def dup_ff_lcm(f, g, K):
h = dup_exquo(dup_mul(f, g, K),
dup_gcd(f, g, K), K)
- return dup_ground_monic(h, K)
+ return dup_monic(h, K)
def dup_lcm(f, g, K):
"""Computes polynomial LCM of `f` and `g` in `K[x]`. """
@@ -1252,6 +1308,7 @@ def dup_lcm(f, g, K):
else:
return dup_rr_lcm(f, g, K)
+@cythonized("u")
def dmp_rr_lcm(f, g, u, K):
"""Computes polynomial LCM over a ring in `K[X]`. """
fc, f = dmp_ground_primitive(f, u, K)
@@ -1264,13 +1321,15 @@ def dmp_rr_lcm(f, g, u, K):
return dmp_mul_ground(h, c, u, K)
+@cythonized("u")
def dmp_ff_lcm(f, g, u, K):
"""Computes polynomial LCM over a field in `K[X]`. """
h = dmp_exquo(dmp_mul(f, g, u, K),
dmp_gcd(f, g, u, K), u, K)
return dmp_ground_monic(h, u, K)
+@cythonized("u")
def dmp_lcm(f, g, u, K):
"""Computes polynomial LCM of `f` and `g` in `K[X]`. """
if not u:
@@ -1298,16 +1357,20 @@ def dup_trunc(f, p, K):
return dup_strip(g)
+@cythonized("u")
def dmp_trunc(f, p, u, K):
"""Reduce `K[X]` polynomial modulo a polynomial `p` in `K[Y]`. """
return dmp_strip([ dmp_rem(c, p, u-1, K) for c in f ], u)
+@cythonized("u,v")
def dmp_ground_trunc(f, p, u, K):
"""Reduce `K[X]` polynomial modulo a constant `p` in `K`. """
if not u:
return dup_trunc(f, p, K)
- else:
- return dmp_strip([ dmp_ground_trunc(c, p, u-1, K) for c in f ], u)
+
+ v = u-1
+
+ return dmp_strip([ dmp_ground_trunc(c, p, v, K) for c in f ], u)
def dup_monic(f, K):
"""Divides all coefficients by `LC(f)` in `K[x]`. """
@@ -1321,6 +1384,7 @@ def dup_monic(f, K):
else:
return dup_quo_ground(f, lc, K)
+@cythonized("u")
def dmp_ground_monic(f, u, K):
"""Divides all coefficients by `LC(f)` in `K[X]`. """
if not u:
@@ -1362,6 +1426,7 @@ def dup_content(f, K):
else:
return dup_rr_content(f, K)
+@cythonized("u,v")
def dmp_content(f, u, K):
"""Returns GCD of multivariate coefficients. """
cont, v = dmp_LC(f, K), u-1
@@ -1380,6 +1445,7 @@ def dmp_content(f, u, K):
else:
return cont
+@cythonized("u,v")
def dmp_rr_ground_content(f, u, K):
"""Returns GCD of coefficients over a ring. """
if not u:
@@ -1396,13 +1462,15 @@ def dmp_rr_ground_content(f, u, K):
return cont
+@cythonized("u")
def dmp_ff_ground_content(f, u, K):
"""Returns GCD of coefficients over a field. """
if not f:
return K.zero
else:
return K.one
+@cythonized("u")
def dmp_ground_content(f, u, K):
"""Returns GCD of coefficients in `K[X]`. """
if not u:
@@ -1433,6 +1501,7 @@ def dup_primitive(f, K):
else:
return dup_rr_primitive(f, K)
+@cythonized("u,v")
def dmp_primitive(f, u, K):
"""Returns multivariate content and a primitive polynomial. """
cont, v = dmp_content(f, u, K), u-1
@@ -1442,6 +1511,7 @@ def dmp_primitive(f, u, K):
else:
return cont, [ dmp_exquo(c, cont, v, K) for c in f ]
+@cythonized("u")
def dmp_rr_ground_primitive(f, u, K):
"""Returns content and a primitive polynomial over a ring. """
cont = dmp_ground_content(f, u, K)
@@ -1451,13 +1521,15 @@ def dmp_rr_ground_primitive(f, u, K):
else:
return cont, dmp_exquo_ground(f, cont, u, K)
+@cythonized("u")
def dmp_ff_ground_primitive(f, u, K):
"""Returns content and a primitive polynomial over a ring. """
if dmp_zero_p(f, u):
return K.zero, f
else:
return K.one, f
+@cythonized("u")
def dmp_ground_primitive(f, u, K):
"""Returns content and a primitive polynomial in `K[x]`. """
if not u:
@@ -1478,6 +1550,7 @@ def dup_sqf_p(f, K):
else: