[PATCH] Fixed subresultant PRS computation and ratint()
Mateusz Paprocki
a side effect, ratint() failed to compute integrals properly. This
patch fixes subresultant algorithm, so that it now returns correct
PRS and resultant (equivalent to result from poly_resultant()).
This fixes ratint() and solves issues #1416, #1417 and #1419.
Signed-off-by: Mateusz Paprocki <mat...@gmail.com>
---
sympy/integrals/rationaltools.py | 10 +++-
sympy/integrals/tests/test_integrals.py | 14 ++----
sympy/integrals/tests/test_rationaltools.py | 28 +++++++++++-
sympy/polys/algorithms.py | 65 +++++++++++++++++++-------
sympy/polys/tests/test_polynomial.py | 28 ++++++++----
sympy/polys/wrappers.py | 18 +++++++-
6 files changed, 121 insertions(+), 42 deletions(-)
diff --git a/sympy/integrals/rationaltools.py b/sympy/integrals/rationaltools.py
index 97657e0..6e3fbc5 100644
--- a/sympy/integrals/rationaltools.py
+++ b/sympy/integrals/rationaltools.py
@@ -40,12 +40,16 @@ def ratint(f, x, **flags):
result, p = poly_div(p, q)
+ result = result.integrate(x).as_basic()
+
+ if p.is_zero:
+ return result
+
g, h = ratint_ratpart(p, q, x)
P, Q = h.as_numer_denom()
q, r = poly_div(P, Q, x)
- result = result.integrate(x).as_basic()
result += g + q.integrate(x).as_basic()
if not r.is_zero:
@@ -150,8 +154,8 @@ def ratint_logpart(f, g, x, t=None):
t = t or Symbol('t', dummy=True)
a, b = g, f - g.diff().mul_term(t)
- R = poly_subresultants(a, b)
- Q = poly_sqf(Poly(R[-1], t))
+ res, R = poly_subresultants(a, b)
+ Q = poly_sqf(Poly(res, t))
R_map, H, i = {}, [], 1
diff --git a/sympy/integrals/tests/test_integrals.py b/sympy/integrals/tests/test_integrals.py
index 301090c..92f4441 100644
--- a/sympy/integrals/tests/test_integrals.py
+++ b/sympy/integrals/tests/test_integrals.py
@@ -146,16 +146,6 @@ def test_issue565():
assert integrate(-Rational(1)/2 * x * sin(n * pi * x/2), [x, -2, 0]) \
== 2*cos(pi*n)/(pi*n)
-
-def test_rational_functions():
- half = Rational(1,2)
- integrate(1/(x**2+x+1), x) == -I*3**half*log(half + x - half*I*3**half)/3 +\
- I*3**half*log(half + x + half*I*3**half)/3
-
- integrate(1/(x**3+1), x) == log(1 + x)/3 - \
- (Rational(1,6) - I*3**half/6)*log(half - x - I*3**half/2) - \
- (Rational(1,6) + I*3**half/6)*log(half - x + I*3**half/2)
-
def test_issue580():
# definite integration of rational functions gives wrong answers
assert NS(Integral(1/(x**2-8*x+17), (x, 2, 4))) == '1.10714871779409'
@@ -182,6 +172,9 @@ def test_issue853():
assert integrate(f, x) == -cos(x)
raises(ValueError, "integrate(f, 2*x)")
+def test_issue1417():
+ assert integrate(2**x - 2*x, x) == 2**x/log(2) - x**2
+
def test_matrices():
M = Matrix(2, 2, lambda i, j: (i+j+1)*sin((i+j+1)*x))
@@ -339,3 +332,4 @@ def test_subs5():
def test_integration_variable():
raises(ValueError, "Integral(exp(-x**2), 3)")
raises(ValueError, "Integral(exp(-x**2), (3, -oo, oo))")
+
diff --git a/sympy/integrals/tests/test_rationaltools.py b/sympy/integrals/tests/test_rationaltools.py
index e212b60..8d33ceb 100644
--- a/sympy/integrals/tests/test_rationaltools.py
+++ b/sympy/integrals/tests/test_rationaltools.py
@@ -1,10 +1,11 @@
-from sympy import symbols, S, I, atan, log
+from sympy import symbols, S, I, atan, log, Poly
from sympy.integrals.rationaltools import ratint, \
ratint_ratpart, ratint_logpart, log_to_atan, log_to_real
x, t = symbols('x t')
+half = S(1)/2
def test_ratint():
f = S(1)
@@ -13,6 +14,11 @@ def test_ratint():
assert ratint(f / g, x) == log(x + 1)
assert ratint((f,g), x) == log(x + 1)
+ f = x**3 - x
+ g = x - 1
+
+ assert ratint(f/g, x) == x**3/3 + x**2/2
+
f = S(1)
g = x**2 + 1
@@ -50,3 +56,23 @@ def test_ratint():
x + S(1)/2*x**2 + S(1)/2*log(2-x+x**2) + (S(9)/7-4*x/7)/(2-x+x**2) + \
13*7**(S(1)/2)*atan(-S(1)/7*7**(S(1)/2) + 2*x*7**(S(1)/2)/7)/49
+ assert ratint(1/(x**2+x+1), x) == \
+ 2*3**(S(1)/2)*atan(3**(S(1)/2)/3 + 2*x*3**(S(1)/2)/3)/3
+
+ assert ratint(1/(x**3+1), x) == \
+ -log(1 - x + x**2)/6 + log(1 + x)/3 + 3**(S(1)/2)*atan(-3**(S(1)/2)/3 + 2*x*3**(S(1)/2)/3)/3
+
+ assert ratint(1/(x**2+x+1), x, real=False) == \
+ -I*3**half*log(half + x - half*I*3**half)/3 + \
+ I*3**half*log(half + x + half*I*3**half)/3
+
+ assert ratint(1/(x**3+1), x, real=False) == log(1 + x)/3 - \
+ (S(1)/6 - I*3**half/6)*log(-half + x + I*3**half/2) - \
+ (S(1)/6 + I*3**half/6)*log(-half + x - I*3**half/2)
+
+def test_ratint_logpart():
+ assert ratint_logpart(x, x**2-9, x, t) == \
+ [(Poly(x**2 - 9, x), Poly(-2*t + 1, t))]
+ assert ratint_logpart(x**2, x**3-5, x, t) == \
+ [(Poly(x**3 - 5, x), Poly(-3*t + 1, t))]
+
diff --git a/sympy/polys/algorithms.py b/sympy/polys/algorithms.py
index b304728..455cf57 100644
--- a/sympy/polys/algorithms.py
+++ b/sympy/polys/algorithms.py
@@ -483,7 +483,7 @@ def poly_gcd(f, g, *symbols):
if f.is_multivariate:
h = poly_div(f*g, poly_lcm(f, g))[0]
else:
- h = poly_subresultants(f, g)[-1]
+ h = poly_subresultants(f, g, res=False)[-1]
if gcd != 1:
return h.mul_term(gcd / h.LC)
@@ -628,7 +628,7 @@ def poly_resultant(f, g, *symbols):
else:
return sign * Poly.cancel(det)
-def poly_subresultants(f, g, *symbols):
+def poly_subresultants(f, g, *symbols, **flags):
"""Computes subresultant PRS of two univariate polynomials.
Polynomial remainder sequence (PRS) is a fundamental tool in
@@ -673,8 +673,8 @@ def poly_subresultants(f, g, *symbols):
if f.is_multivariate:
raise MultivariatePolyError(f)
-
- symbols, flags = f.symbols, f.flags
+ else:
+ symbols = f.symbols
n, m = f.degree, g.degree
@@ -682,36 +682,65 @@ def poly_subresultants(f, g, *symbols):
f, g = g, f
n, m = m, n
- prs = [f, g]
+ R = [f, g]
d = n - m
- b = (-1)**(d + 1)
-
- h = poly_pdiv(f, g)[1]
- h = h.mul_term(b)
+ b = S(-1)**(d + 1)
+ c = S(-1)
- k = h.degree
+ B, D = [b], [d]
- c = S.NegativeOne
+ h = poly_prem(f, g)
+ h = h.mul_term(b)
while not h.is_zero:
- prs.append(h)
+ k = h.degree
+ R.append(h)
- coeff = g.LC
+ lc = g.LC
- c = (-coeff)**d / c**(d-1)
+ C = (-lc)**d / c**(d-1)
+ c = Poly.cancel(C)
- b = -coeff * c**(m-k)
+ b = -lc * c**(m-k)
f, g, m, d = g, h, k, m-k
- h = poly_pdiv(f, g)[1]
+ B.append(b)
+ D.append(d)
+
+ h = poly_prem(f, g)
h = h.div_term(b)
- k = h.degree
+ if not flags.get('res', True):
+ return R
+
+ if R[-1].degree > 0:
+ return (Poly((), *symbols), R)
+ if R[-2].is_one:
+ return (R[-1], R)
+
+ s, c, i = 1, S(1), 1
+
+ for b, d in zip(B, D)[:-1]:
+ u = R[i-1].degree
+ v = R[i ].degree
+ w = R[i+1].degree
+
+ if u % 2 and v % 2:
+ s = -s
+
+ lc = R[i].LC
+
+ C = c*(b/lc**(1 + d))**v * lc**(u - w)
+ c = Poly.cancel(C)
+
+ i += 1
+
+ j = R[-2].degree
- return prs
+ return (R[-1]**j*s*c, R)
def poly_sqf(f, *symbols):
"""Compute square-free decomposition of an univariate polynomial.
diff --git a/sympy/polys/tests/test_polynomial.py b/sympy/polys/tests/test_polynomial.py
index 4dbc9cf..7989bb6 100644
--- a/sympy/polys/tests/test_polynomial.py
+++ b/sympy/polys/tests/test_polynomial.py
@@ -615,14 +615,24 @@ def test_poly_subresultants():
g = 3*x**6+5*x**4-4*x**2-9*x+21
assert poly_subresultants(f, g, x) == \
- [Poly(f, x), Poly(g, x),
- Poly(15*x**4 - 3*x**2 + 9, x),
- Poly(65*x**2 + 125*x - 245, x),
- Poly(9326*x - 12300, x),
- Poly(260708, x)]
+ (Poly(260708, x), [Poly(f, x), Poly(g, x),
+ Poly(15*x**4 - 3*x**2 + 9, x),
+ Poly(65*x**2 + 125*x - 245, x),
+ Poly(9326*x - 12300, x),
+ Poly(260708, x)])
assert poly_subresultants((x-1)**2, x**2-1, x) == \
- [Poly((x-1)**2, x), Poly(x**2-1, x), Poly(2*x - 2, x)]
+ (Poly(0, x), [Poly((x-1)**2, x),
+ Poly(x**2-1, x),
+ Poly(2*x - 2, x)])
+
+ f = Poly(-x**3 + 5, x)
+ g = Poly((1 + 3*t)*x**2, x)
+
+ assert poly_subresultants(f, g) == \
+ (Poly(25 + 225*t + 675*t**2 + 675*t**3, x), [Poly(-x**3 + 5, x),
+ Poly((1 + 3*t)*x**2, x),
+ Poly(5 + 30*t + 45*t**2, x)])
def test_poly_groebner():
assert poly_groebner(0, x) == [Poly((), x)]
@@ -773,10 +783,10 @@ def test_squarefree():
A = Poly(x**4 - 3*x**2 + 6, x)
D = Poly(x**6 - 5*x**4 + 5*x**2 + 4, x)
- f, g = D, A - D.diff(x)*t
+ f, g = D, A - D.diff(x).mul_term(t)
- R = poly_subresultants(f, g)
- S = poly_sqf(Poly(R[-1], t))
+ res, R = poly_subresultants(f, g)
+ S = poly_sqf(Poly(res, t))
assert S == [Poly(45796, t), Poly(1, t), Poly(4*t**2 + 1, t)]
diff --git a/sympy/polys/wrappers.py b/sympy/polys/wrappers.py
index 44dc1b8..2829161 100644
--- a/sympy/polys/wrappers.py
+++ b/sympy/polys/wrappers.py
@@ -38,7 +38,6 @@ def _map_basic(f, n, *args, **kwargs):
'gcd' : 2,
'gcdex' : 2,
'half_gcdex' : 2,
- 'subresultants' : 2,
'resultant' : 2,
'sqf' : 1,
'decompose' : 1,
@@ -69,3 +68,20 @@ def div(*args, **kwargs):
div.__doc__ = poly_div.__doc__
+def subresultants(*args, **kwargs):
+ result = poly_subresultants(*_conv_args(2, args), **kwargs)
+
+ if type(result) is tuple:
+ res, R = result
+ else:
+ res, R = None, result
+
+ R = [ r.as_basic() for r in R ]
+
+ if res is None:
+ return R
+ else:
+ return res.as_basic(), R
+
+subresultants.__doc__ = poly_subresultants.__doc__
+
--
1.6.2.3
Ondrej Certik
Please push it in.
Ondrej
Jochen Voß
I left a few comments about your patch in the "Fix traceback when
calculating simple integrals" thread. Sorry for misplacing them.
All the best,
Jochen
--
http://seehuhn.de/