The following two patches address in one patch several issues that
I've raised in various issues. I've tried to keep the patch to items
that are absolutely essential to change to lighten the reviewer's
burden. Addressed are:
the quartic (and cubic) root finding routines and changes to roots to
be compatible with the changes;
changes to avoid the infinite recursion that would be generated by (x*
(-x-x**3)).could_extract_minus_sign()
fixes to -base (int or rat) ** -exp (int or rat)
changes to how complex multiplication handles uncertainties
I left out most string modifications that were put in other patches
that are attached to the issues.
Tests have also been added and modified as necessary.
/c
========= patch that adds one more test
>From 6d663d62e3bad7bbaaae47b2233f57618b220a74 Mon Sep 17 00:00:00 2001
From: Chris Smith <smi...@gmail.com>
Date: Sun, 2 Aug 2009 15:46:51 +0545
Subject: [PATCH 1/1] test for infinite recursive minus extraction
---
sympy/core/tests/test_basic.py | 1 +
1 files changed, 1 insertions(+), 0 deletions(-)
diff --git a/sympy/core/tests/test_basic.py b/sympy/core/tests/
test_basic.py
index 7c51702..e0226ac 100644
--- a/sympy/core/tests/test_basic.py
+++ b/sympy/core/tests/test_basic.py
@@ -593,6 +593,7 @@ def test_extractions():
assert (-(x+x*y)/y).could_extract_minus_sign() == True
assert ((x+x*y)/(-y)).could_extract_minus_sign() == True
assert ((x+x*y)/y).could_extract_minus_sign() == False
+ assert (x*(-x-x**3)).could_extract_minus_sign() == True
def test_coeff():
assert (3+2*x+4*x**2).coeff(1) == None
--
1.6.3.2.1299.gee46c
======== patch that adds everything else
>From 432fe52a266ffb5dede6656bb44f51982efacb45 Mon Sep 17 00:00:00 2001
From: Chris Smith <smi...@gmail.com>
Date: Sat, 1 Aug 2009 19:00:21 +0545
Subject: [PATCH 1/1] re,im_acc switch made in evalf_mul to handle
problems with
complex evalf losing precision, but that didn't fix everything.
The complete change that gets closer (hopefully) to the right
things along with other critical updates to documented problems
are a part of this patch:
replaced the quartic solver with a new one that solves what
the other didn't
made the cubic routine complete in its solutions (it was
not set up to give solutions for cases that should
be picked up by roots()
made roots take a list of coefficients to represent the
coefficients of a univariate polynomial with the last
one representing the constant term. This is needed for the
updates to the quartic and cubic routines.
Fixed the infinite recurssion when trying to extract a minus
sign by changing how basic.could_extract_minus_sign() works
for the is_Mul case. In the is_Add case, we are giving preference
to an expression that has more positives than negatives. The code
for the is_Mul case was just breaking the expr into num and den
and then *if the den was not zero" (which only happens for oo)
the numerator was checked to see if a sign could be extracted. The
peoblem is, however, that this led to infinete recursion if the
expr didn't have a denominator since the numerator just kept getting
recycled into the routine. I think what we should do is check the
num and den and give preference to the expression that has an odd
number of negatives that could be extracted. This would be consistent
with an expression having a negative in front of it once all other
terms were in standard/canonical form.
The exponentiation of negative Int and Rat by negative powers
(Both Int and Rat) has been fixed.
---
sympy/concrete/tests/test_gosper.py | 2 +-
sympy/core/__init__.py | 4 +-
sympy/core/basic.py | 8 +-
sympy/core/evalf.py | 53 +++++------
sympy/core/numbers.py | 15 ++-
sympy/core/tests/test_basic.py | 9 ++
sympy/core/tests/test_evalf.py | 7 +-
sympy/core/tests/test_functions.py | 2 +-
sympy/core/tests/test_multidimensional.py | 2 +-
sympy/core/tests/test_numbers.py | 16 ++--
sympy/functions/combinatorial/__init__.py | 2 +-
sympy/mpmath/optimization.py | 14 ++--
sympy/polys/polynomial.py | 22 ++++-
sympy/polys/rootfinding.py | 145 +++++++++++++++++
+-----------
sympy/polys/tests/test_polynomial.py | 92 +++++++++++--------
sympy/simplify/tests/test_simplify.py | 2 +
16 files changed, 243 insertions(+), 152 deletions(-)
diff --git a/sympy/concrete/tests/test_gosper.py b/sympy/concrete/
tests/test_gosper.py
index bd07751..5887fb9 100644
--- a/sympy/concrete/tests/test_gosper.py
+++ b/sympy/concrete/tests/test_gosper.py
@@ -3,4 +3,4 @@ def test_normal():
pass
def test_gosper():
- pass
\ No newline at end of file
+ pass
diff --git a/sympy/core/__init__.py b/sympy/core/__init__.py
index 516df7e..635fc24 100644
--- a/sympy/core/__init__.py
+++ b/sympy/core/__init__.py
@@ -3,8 +3,8 @@
from basic import Basic, S, C, sympify
from symbol import Symbol, Wild, symbols, var
-from numbers import Number, Real, Rational, Integer, igcd, ilcm,
RealNumber, \
- seterr
+from numbers import Number, Real, Rational, Integer, NumberSymbol, \
+ RealNumber, igcd, ilcm, seterr
from power import Pow, integer_nthroot
from mul import Mul
from add import Add
diff --git a/sympy/core/basic.py b/sympy/core/basic.py
index b20ed06..d336623 100644
--- a/sympy/core/basic.py
+++ b/sympy/core/basic.py
@@ -2014,9 +2014,13 @@ def could_extract_minus_sign(self):
elif positive_args < negative_args:
return True
elif self.is_Mul:
+ # We choose the one with an odd number of minus signs
num, den = self.as_numer_denom()
- if den != 0:
- return num.could_extract_minus_sign()
+ args = (list(num.args) if num.is_Mul else [num]) + \
+ (list(den.args) if den.is_Mul else [den])
+ arg_signs = [arg.could_extract_minus_sign() for arg
in args]
+ negative_args = filter(None, arg_signs)
+ return len(negative_args) % 2 == 1
# As a last resort, we choose the one with greater hash
return hash(self) < hash(negative_self)
diff --git a/sympy/core/evalf.py b/sympy/core/evalf.py
index fd72fca..7846484 100644
--- a/sympy/core/evalf.py
+++ b/sympy/core/evalf.py
@@ -312,18 +312,6 @@ def evalf_add(v, prec, options):
finally:
options['maxprec'] = oldmaxprec
-# Helper for complex multiplication
-# XXX: should be able to multiply directly, and use complex_accuracy
-# to obtain the final accuracy
-def cmul((a, aacc), (b, bacc), (c, cacc), (d, dacc), prec,
target_prec):
- A, Aacc = mpf_mul(a,c,prec), min(aacc, cacc)
- B, Bacc = mpf_mul(mpf_neg(b),d,prec), min(bacc, dacc)
- C, Cacc = mpf_mul(a,d,prec), min(aacc, dacc)
- D, Dacc = mpf_mul(b,c,prec), min(bacc, cacc)
- re, re_accuracy = add_terms([(A, Aacc), (B, Bacc)], prec,
target_prec)
- im, im_accuracy = add_terms([(C, Cacc), (D, Cacc)], prec,
target_prec)
- return re, im, re_accuracy, im_accuracy
-
def evalf_mul(v, prec, options):
args = v.args
# With guard digits, multiplication in the real case does not
destroy
@@ -343,16 +331,15 @@ def evalf_mul(v, prec, options):
# direction tells us that the result should be multiplied by
# i**direction
for arg in args:
- re, im, a, aim = evalf(arg, prec, options)
+ re, im, re_acc, im_acc = evalf(arg, prec, options)
if re and im:
- complex_factors.append((re, im, a, aim))
+ complex_factors.append((re, im, re_acc, im_acc))
continue
elif re:
- s, m, e, b = re
+ (s, m, e, b), w_acc = re, re_acc
elif im:
- a = aim
+ (s, m, e, b), w_acc = im, im_acc
direction += 1
- s, m, e, b = im
else:
return None, None, None, None
direction += 2*s
@@ -362,28 +349,32 @@ def evalf_mul(v, prec, options):
if bc > 3*prec:
man >>= prec
exp += prec
- acc = min(acc, a)
+ acc = min(acc, w_acc)
sign = (direction & 2) >> 1
v = normalize(sign, man, exp, bitcount(man), prec, round_nearest)
if complex_factors:
- # Multiply first complex number by the existing real scalar
- re, im, re_acc, im_acc = complex_factors[0]
- re = mpf_mul(re, v, prec)
- im = mpf_mul(im, v, prec)
- re_acc = min(re_acc, acc)
- im_acc = min(im_acc, acc)
- # Multiply consecutive complex factors
- complex_factors = complex_factors[1:]
+ # make existing real scalar look like an imaginary and
+ # multiply by the remaining complex numbers
+ re, im = v, (0, MP_BASE(0), 0, 0)
for wre, wim, wre_acc, wim_acc in complex_factors:
- re, im, re_acc, im_acc = cmul((re, re_acc), (im,im_acc),
- (wre,wre_acc), (wim,wim_acc), prec, target_prec)
+ # acc is the overall accuracy of the product; we aren't
+ # computing exact accuracies of the product.
+ acc = min(acc,
+ complex_accuracy((wre, wim, wre_acc, wim_acc)))
+ A = mpf_mul(re, wre, prec)
+ B = mpf_mul(mpf_neg(im), wim, prec)
+ C = mpf_mul(re, wim, prec)
+ D = mpf_mul(im, wre, prec)
+ re, xre_acc = add_terms([(A, acc), (B, acc)], prec,
target_prec)
+ im, xim_acc = add_terms([(C, acc), (D, acc)], prec,
target_prec)
+
if options.get('verbose'):
- print "MUL: obtained accuracy", re_acc, im_acc,
"expected", target_prec
+ print "MUL: wanted", target_prec, "accurate bits, got",
acc
# multiply by i
if direction & 1:
- return mpf_neg(im), re, re_acc, im_acc
+ return mpf_neg(im), re, acc, acc
else:
- return re, im, re_acc, im_acc
+ return re, im, acc, acc
else:
# multiply by i
if direction & 1:
diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py
index 11a2b72..8b1a9fc 100644
--- a/sympy/core/numbers.py
+++ b/sympy/core/numbers.py
@@ -583,7 +583,10 @@ def _eval_power(b, e):
if (ne is S.One):
return Rational(b.q, b.p)
if b < 0:
- return -(S.NegativeOne) ** ((e.p % e.q) / S(e.q))
* Rational(b.q, -b.p) ** ne
+ if e.q != 1:
+ return -(S.NegativeOne) ** ((e.p % e.q) / S
(e.q)) * Rational(b.q, -b.p) ** ne
+ else:
+ return S.NegativeOne ** ne * Rational(b.q, -
b.p) ** ne
else:
return Rational(b.q, b.p) ** ne
if (e is S.Infinity):
@@ -650,7 +653,7 @@ def __ne__(self, other):
if other.is_comparable and not isinstance(other, Rational):
other = other.evalf()
if isinstance(other, Number):
if isinstance(other, Real):
- return bool(not mlib.feq(self._as_mpf_val
(other._prec), other._mpf_))
+ return bool(not mlib.mpf_eq(self._as_mpf_val
(other._prec), other._mpf_))
return bool(self.p!=other.p or self.q!=other.q)
return True # Rational != non-Number
@@ -944,10 +947,14 @@ def _eval_power(base, exp):
return S.ImaginaryUnit * Pow(-base, exp)
if exp < 0:
# invert base and change sign on exponent
+ ne = -exp
if base < 0:
- return -(S.NegativeOne) ** ((exp.p % exp.q) / S
(exp.q)) * Rational(1, -base) ** (-exp)
+ if exp.q != 1:
+ return -(S.NegativeOne) ** ((exp.p % exp.q) / S
(exp.q)) * Rational(1, -base) ** ne
+ else:
+ return (S.NegativeOne) ** ne * Rational(1, -base)
** ne
else:
- return Rational(1, base.p) ** (-exp)
+ return Rational(1, base.p) ** ne
# see if base is a perfect root, sqrt(4) --> 2
x, xexact = integer_nthroot(abs(base.p), exp.q)
if xexact:
diff --git a/sympy/core/tests/test_basic.py b/sympy/core/tests/
test_basic.py
index 7c51702..6b2a1c7 100644
--- a/sympy/core/tests/test_basic.py
+++ b/sympy/core/tests/test_basic.py
@@ -361,6 +361,12 @@ def test_noncommutative_expand_issue658():
assert (A*(A+B)*B).expand() == A**2*B + A*B**2
assert (A*(A+B+C)*B).expand() == A**2*B + A*B**2 + A*C*B
+def test_as_numer_denom():
+ assert oo.as_numer_denom() == (1, 0)
+ assert (1/x).as_numer_denom() == (1, x)
+ assert x.as_numer_denom() == (x, 1)
+ assert (x/-y).as_numer_denom() == (-x, y)
+
def test_as_independent():
assert (2*x*sin(x)+y+x).as_independent(x) == (y, x + 2*x*sin(x))
assert (2*x*sin(x)+y+x).as_independent(y) == (x + 2*x*sin(x), y)
@@ -593,6 +599,9 @@ def test_extractions():
assert (-(x+x*y)/y).could_extract_minus_sign() == True
assert ((x+x*y)/(-y)).could_extract_minus_sign() == True
assert ((x+x*y)/y).could_extract_minus_sign() == False
+ assert (x*(-x - x**2)).could_extract_minus_sign() == True
+ assert ((-x-y)/(x+y)).could_extract_minus_sign() == True #is_Mul
odd case
+ assert ((-x-y)/(x-y)).could_extract_minus_sign() == False #is_Mul
even case
def test_coeff():
assert (3+2*x+4*x**2).coeff(1) == None
diff --git a/sympy/core/tests/test_evalf.py b/sympy/core/tests/
test_evalf.py
index 4782f45..385b3d9 100644
--- a/sympy/core/tests/test_evalf.py
+++ b/sympy/core/tests/test_evalf.py
@@ -87,8 +87,8 @@ def test_evalf_complex_cancellation():
# XXX: the number of returned mantissa digits in the real part
could
# change with the implementation. What matters is that the
returned digits are
# correct.
- assert NS((A+B*I)*(C+D*I),6) in ('6.45e-6 + 0.892529*I', '6.4e-6
+ 0.892529*I')
- assert NS((A+B*I)*(C+D*I),10) == '6.447100e-6 + 0.8925286452*I'
+ assert NS((A+B*I)*(C+D*I),6) == '6.44862e-6 + 0.892529*I'
+ assert NS((A+B*I)*(C+D*I),10) == '6.447099821e-6 +
0.8925286452*I'
assert NS((A+B*I)*(C+D*I) - F*I, 5) in ('6.4471e-6 - .0e-15*I',
'6.4471e-6 + .0e-15*I')
def test_evalf_logs():
@@ -137,6 +137,9 @@ def test_evalf_bugs():
assert NS('log(2)',10) == '0.6931471806'
assert NS('(sin(x)-x)/x**3', 15, subs={x:'1/10**50'}) ==
'-0.166666666666667'
assert NS(sin(1)+Rational(1,10**100)*I,15) == '0.841470984807897
+ 1.00000000000000e-100*I'
+ assert NS((1+I)**2*I,6) == '-2.00000 + 2.32831e-10*I'
+ d={n: (-1)**Rational(6,7), y: (-1)**Rational(4,7), x: (-1)
**Rational(2,7)}
+ assert NS((x*(1+y*(1 + n))).subs(d).evalf(),6) == '0.346011 +
0.433884*I'
assert x.evalf() == x
def test_evalf_integer_parts():
diff --git a/sympy/core/tests/test_functions.py b/sympy/core/tests/
test_functions.py
index a1cdf10..3e99a55 100644
--- a/sympy/core/tests/test_functions.py
+++ b/sympy/core/tests/test_functions.py
@@ -295,4 +295,4 @@ class MyFunc(Function):
@classmethod
def eval(cls, *args):
return (0,0,0)
- assert MyFunc(0) == (0,0,0)
\ No newline at end of file
+ assert MyFunc(0) == (0,0,0)
diff --git a/sympy/core/tests/test_multidimensional.py b/sympy/core/
tests/test_multidimensional.py
index 6e26602..fb06a56 100644
--- a/sympy/core/tests/test_multidimensional.py
+++ b/sympy/core/tests/test_multidimensional.py
@@ -38,4 +38,4 @@ def test_diffmulti():
def test_expandmulti():
x = Symbol("x")
y = Symbol("y")
- assert expand([(x+y)**2, (x-1)**2])==[expand((x+y)**2), expand
((x-1)**2)]
\ No newline at end of file
+ assert expand([(x+y)**2, (x-1)**2])==[expand((x+y)**2), expand
((x-1)**2)]
diff --git a/sympy/core/tests/test_numbers.py b/sympy/core/tests/
test_numbers.py
index 876d782..cffe2ca 100644
--- a/sympy/core/tests/test_numbers.py
+++ b/sympy/core/tests/test_numbers.py
@@ -97,7 +97,7 @@ def test_Rational_cmp():
assert not n3<n1
assert not (Rational(-1) > 0)
assert Rational(-1) < 0
-
+ assert S(1) != 2**(S(1)/2)
def test_Real():
def eq(a, b):
@@ -221,6 +221,7 @@ def test_powers_Integer():
assert S(-9) ** Rational(3, 2) == -27*I
assert S(27) ** Rational(2, 3) == 9
assert S(-27) ** Rational(2, 3) == 9 * (S(-1) ** Rational(2, 3))
+ assert (-2) ** Rational(-2, 1) == Rational(1, 4)
# not exact roots
assert (-3) ** (S(1)/2) == sqrt(-3)
@@ -232,6 +233,8 @@ def test_powers_Integer():
assert (2) ** (S(-3)/2) == sqrt(2) / 4
assert (81) ** (S(2)/3) == 9 * (S(3) ** (S(2)/3))
assert (-81) ** (S(2)/3) == 9 * (S(-3) ** (S(2)/3))
+ assert (-3) ** Rational(-7, 3) == -(-3) ** Rational(2, 3) / 27
+ assert (-3) ** Rational(-2, 3) == -(-3) ** Rational(1, 3) / 3
# join roots
@@ -247,13 +250,9 @@ def test_powers_Integer():
assert (2**64+1) ** Rational(4, 3)
assert (2**64+1) ** Rational(17,25)
- # negative rational power and negative base
- assert (-3) ** Rational(-7, 3) == -(-3) ** Rational(2, 3) / 27
- assert (-3) ** Rational(-2, 3) == -(-3) ** (S(1) / 3) / 3
-
def test_powers_Rational():
"""Test Rational._eval_power"""
- # check inifinity
+ # check infinity
assert Rational(1,2) ** S.Infinity == 0
assert Rational(3,2) ** S.Infinity == S.Infinity
assert Rational(-1,2) ** S.Infinity == 0
@@ -282,12 +281,15 @@ def test_powers_Rational():
assert Rational(1,2) ** Rational(1,2) == sqrt(2) / 2
assert Rational(-4,7) ** Rational(1,2) == I * Rational(4,7) **
Rational(1,2)
- # negative rational power and negative base
+ # negative rational power and negative rational base
assert Rational(-3, 2)**Rational(-7, 3) == \
-4 * (-3) ** Rational(2, 3)*2 ** Rational(1, 3)/27
assert Rational(-3, 2)**Rational(-2, 3) == \
-(-3) ** (S(1) / 3) * 2 ** (S(2) / 3) / 3
+ # negative integer power and negative rational base
+ assert Rational(-2, 3)**Rational(-2, 1) == Rational(9, 4)
+
def test_abs1():
assert Rational(1,6) != Rational(-1,6)
assert abs(Rational(1,6)) == abs(Rational(-1,6))
diff --git a/sympy/functions/combinatorial/__init__.py b/sympy/
functions/combinatorial/__init__.py
index 8271fc2..78c96d7 100644
--- a/sympy/functions/combinatorial/__init__.py
+++ b/sympy/functions/combinatorial/__init__.py
@@ -1,3 +1,3 @@
import factorials
-import numbers
\ No newline at end of file
+import numbers
diff --git a/sympy/mpmath/optimization.py b/sympy/mpmath/
optimization.py
index 1834ab2..ef6b2ec 100644
--- a/sympy/mpmath/optimization.py
+++ b/sympy/mpmath/optimization.py
@@ -32,7 +32,7 @@ def __init__(self, f, x0, **kwargs):
if len(x0) == 1:
self.x0 = x0[0]
else:
- raise ValueError('expected 1 starting point, got %i' * len
(x0))
+ raise ValueError('expected 1 starting point, got %i' % len
(x0))
self.f = f
if not 'df' in kwargs:
def df(x):
@@ -76,7 +76,7 @@ def __init__(self, f, x0, **kwargs):
self.x0 = x0[0]
self.x1 = x0[1]
else:
- raise ValueError('expected 1 or 2 starting points, got
%i' * len(x0))
+ raise ValueError('expected 1 or 2 starting points, got
%i' % len(x0))
self.f = f
def __iter__(self):
@@ -117,7 +117,7 @@ class MNewton:
def __init__(self, f, x0, **kwargs):
if not len(x0) == 1:
- raise ValueError('expected 1 starting point, got %i' * len
(x0))
+ raise ValueError('expected 1 starting point, got %i' % len
(x0))
self.x0 = x0[0]
self.f = f
if not 'df' in kwargs:
@@ -299,7 +299,7 @@ class Bisection:
def __init__(self, f, x0, **kwargs):
if len(x0) != 2:
- raise ValueError('expected interval of 2 points, got %i'
* len(x0))
+ raise ValueError('expected interval of 2 points, got %i'
% len(x0))
self.f = f
self.a = x0[0]
self.b = x0[1]
@@ -375,7 +375,7 @@ class Illinois:
def __init__(self, f, x0, **kwargs):
if len(x0) != 2:
- raise ValueError('expected interval of 2 points, got %i'
* len(x0))
+ raise ValueError('expected interval of 2 points, got %i'
% len(x0))
self.a = x0[0]
self.b = x0[1]
self.f = f
@@ -467,7 +467,7 @@ class Ridder:
def __init__(self, f, x0, **kwargs):
self.f = f
if len(x0) != 2:
- raise ValueError('expected interval of 2 points, got %i'
* len(x0))
+ raise ValueError('expected interval of 2 points, got %i'
% len(x0))
self.x1 = x0[0]
self.x2 = x0[1]
self.verbose = kwargs['verbose']
@@ -510,7 +510,7 @@ class ANewton:
def __init__(self, f, x0, **kwargs):
if not len(x0) == 1:
- raise ValueError('expected 1 starting point, got %i' * len
(x0))
+ raise ValueError('expected 1 starting point, got %i' % len
(x0))
self.x0 = x0[0]
self.f = f
if not 'df' in kwargs:
diff --git a/sympy/polys/polynomial.py b/sympy/polys/polynomial.py
index 3648c30..7fc401b 100644
--- a/sympy/polys/polynomial.py
+++ b/sympy/polys/polynomial.py
@@ -88,12 +88,30 @@ class Poly(Basic):
a collection of coefficients and monomials. This can be done
in at least three different ways:
- [1] [(c_1, M_1), (c_2, M_2), ..., (c_1, M_1)]
+ [1] [(c_1, M_1), (c_2, M_2), ..., (c_1, M_n)]
- [2] (c_1, c2, ..., c_n), (M_1, M_2, ..., M_n)
+ # univariate monomials don't have to be tuples
+ >>> Poly([(1, 2), (2, 1), (3, 0)], x)
+ Poly(x**2 + 2*x + 3, x)
+ >>> Poly([(3, (2,1))], x, y)
+ Poly(3*x**2*y, x, y)
+ >>> Poly(zip((1,2,3), reversed(range(3))), x)
+ Poly(x**2 + 2*x + 3, x)
+
+ [2] ((c_1, c_2, ..., c_n), (M_1, M_2, ..., M_n))
+
+ >>> Poly(((S(3), S(2)),((3,),(4,))), x)
+ Poly(3*x**3 + 2*x**4, x)
+ >>> Poly(((S(3), ),((2, ),)), x)
+ Poly(3*x**2, x)
+ >>> Poly(( sympify([3,2,1]), ((i,) for i in range
(3) )), x)
+ Poly(3 + 2*x + x**2, x)
[3] { M_1 : c_1, M_2 : c_2, ..., M_n : c_n }
+ >>> Poly( {(1, 2): S(3), (4, 5): S(6)} ,x ,y)
+ Poly(6*x**4*y**5 + 3*x*y**2, x, y)
+
Although all three representation look similar, they are
designed for different tasks and have specific properties:
diff --git a/sympy/polys/rootfinding.py b/sympy/polys/rootfinding.py
index 26cffd2..05a9d6b 100644
--- a/sympy/polys/rootfinding.py
+++ b/sympy/polys/rootfinding.py
@@ -1,4 +1,5 @@
+from sympy.core.symbol import Symbol
from sympy.core.add import Add
from sympy.core.mul import Mul
from sympy.core.basic import Basic, S
@@ -31,84 +32,114 @@ def roots_quadratic(f):
(-b - d) / (2*a),
]
- from sympy.simplify import simplify
-
- return [ simplify(r) for r in roots ]
+ return roots
def roots_cubic(f):
- """Returns a list of roots of a cubic polynomial."""
- if isinstance(f, (tuple, list)):
- a, b, c, d = f
- else:
- a, b, c, d = f.iter_all_coeffs()
+ """Returns a list of roots of a cubic polynomial."""
+ one, a, b, c = f.as_monic().iter_all_coeffs()
- b, c, d = b/a, c/a, d/a
+ if c is S.Zero:
+ x1, x2 = roots((1,a,b), multiple = True)
+ return [x1, S.Zero, x2]
- p = c - b**2 / 3
- q = d + (2*b**3 - 9*b*c) / 27
+ aon3 = a/3
+ p = b - 3*aon3**2
+ q = c - b*aon3 + 2*aon3**3
+ pon3 = p/3
if p is S.Zero:
if q is S.Zero:
- return [-b/3] * 3
+ return [-aon3] * 3
else:
u1 = q**Rational(1, 3)
+ elif q is S.Zero:
+ y1, y2 = roots((1, 0, p), multiple=True)
+ return [tmp - aon3 for tmp in [y1, S.Zero, y2]]
else:
- u1 = (q/2 + (q**2/4 + p**3/27)**S.Half)**Rational(1, 3)
+ u1 = (q/2 + (q**2/4 + pon3**3)**S.Half)**Rational(1, 3)
coeff = S.ImaginaryUnit*3**S.Half / 2
u2 = u1*(-S.Half + coeff)
u3 = u1*(-S.Half - coeff)
- roots = [
- (p/(3*u1) - u1 - b/3),
- (p/(3*u2) - u2 - b/3),
- (p/(3*u3) - u3 - b/3),
+ soln = [
+ -u1 + pon3/u1 - aon3,
+ -u2 + pon3/u2 - aon3,
+ -u3 + pon3/u3 - aon3
]
- return [ r.expand() for r in roots ]
+ return soln
def roots_quartic(f):
"""Returns a list of roots of a quartic polynomial.
- References:
-
http://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormula.html
- """
- a, b, c, d, e = f.iter_all_coeffs()
-
- # normalize
- t = a
- a = b/t
- b = c/t
- c = d/t
- d = e/t
-
- p = -2*b
- q = b**2 + a*c - 4*d
- r = c**2 + a**2*d - a*b*c
-
- r = roots_cubic((1, p, q, r))
+ There are many references for solving quartic expressions
available [1-6].
+ This reviewer has found that many of them require one to select
from among
+ 2 or more possible sets of solutions and that some solutions work
when one
+ is searching for real roots but don't work when searching for
complex roots
+ (though this is not always stated clearly). The following routine
has been
+ tested and found to be correct for 0, 2 or 4 complex roots.
- u = r[1] + r[2] - r[0]
- v = (u**2 - 16*d)**S.Half
- w = (a**2 - 4*r[0])**S.Half
+ Example::
+ >>> from sympy import *
+ >>> x = var('x')
+ >>> r = rootfinding.roots_quartic(Poly(x**4 -6*x**3 +17*x**2
-26*x +20, x))
+ >>> [tmp.evalf(n=2) for tmp in r] # 4 complex roots
+ [1.0 + 1.7*I, 1.0 - 1.7*I, 2.0 + I, 2.0 - 1.0*I] #1+-I*sqrt(3),
2+-I
- A = (u + v) / 4
- B = (u - v) / 4
- C = (-a + w) / 2
- D = (-a - w) / 2
+ References:
+ [1] http://mathforum.org/dr.math/faq/faq.cubic.equations.html
+ [2]
http://en.wikipedia.org/wiki/Quartic_formula#Solving_a_quartic_equation
+ [3]
http://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormula.html
+ [4] Maxima 0.7.3a solution to z**4+e*z**2+f*z+g=0 where
+ x=z-a/4 and x^4+a^x^3+b^x^2+c*x+d=0
+ [5] http://staff.bath.ac.uk/masjhd/JHD-CA.pdf
+ [6] http://www.albmath.org/files/Math_5713.pdf
- E = (C**2 - 4*A)**S.Half
- F = (D**2 - 4*B)**S.Half
+ """
- roots = [
- (C + E) / 2,
- (C - E) / 2,
- (D + F) / 2,
- (D - F) / 2,
- ]
+ #normalized coefficients
+ one, a, b, c, d = f.as_monic().iter_all_coeffs()
- return [ r.expand() for r in roots ]
+ if d is S.Zero:
+ return [S.Zero] + roots((1, a, b, c), multiple = True)
+ else:
+ a2 = a ** 2
+ e = b - 3 * a2 / 8
+ f = c + a * (a2 / 8 - b / 2)
+ g = d - a * (a * (3 * a2 / 256 - b / 16) + c / 4)
+ aon4 = a / 4
+
+ if f is S.Zero:
+ y1, y2 = [tmp ** S.Half for tmp in
+ roots((1, e, g), multiple = True)]
+ return [tmp - aon4 for tmp in [-y1, -y2, y1, y2]]
+ elif g is S.Zero:
+ y = [S.Zero] + roots((1, 0, e, f), multiple = True)
+ return [tmp - aon4 for tmp in y]
+ else:
+ e2, f2 = e**2, f**2
+ a0 = (3 ** (-S(3) / 2) *\
+ (16 * g * (8 * g * (-2 * g + e2) - e * (9 * f2 + e
* e2)) +\
+ f2 * (27 * f2 + 4 * e * e2)) ** S.Half) / 2
+ a1 = (a0 + (2 * e * (-36 * g + e2) + 27 * f2) / 54) ** (S
(1) / 3)
+ a2 = e * 6 * a1
+ a3 = 9 * a1 * a1 + 12 * g + e2
+ a4 = ((a3 - a2) / a1) ** S.Half
+ a5 = (a3 + 2 * a2)
+ a6 = a1 * f * 54
+ a54 = a5 * a4
+ a7 = (-(a54 - a6) / a1) ** S.Half / (6 * a4 ** S.Half)
+ a8 = (-(a54 + a6) / a1 / a4) ** S.Half / 6
+ a4on6 = a4 / 6
+ soln = [-a7 - a4on6 - aon4,
+ a7 - a4on6 - aon4,
+ -a8 + a4on6 - aon4,
+ a8 + a4on6 - aon4
+ ]
+
+ return soln
def roots_binomial(f):
"""Returns a list of roots of a binomial polynomial."""
@@ -183,12 +214,20 @@ def roots(f, *symbols, **flags):
>>> roots(x**2 - 1, x)
{1: 1, -1: 1}
+ >>> roots((1, 0, -1))
+ {1: 1, -1: 1}
+
>>> roots(x**2 - y, x)
{y**(1/2): 1, -y**(1/2): 1}
"""
if not isinstance(f, Poly):
- f = Poly(f, *symbols)
+ if isinstance(f, (list, tuple)):
+ c = sympify(f)
+ m = ((i,) for i in reversed(range(len(c))))
+ f = Poly((c, m), Symbol('x'))
+ else:
+ f = Poly(f, *symbols)
elif symbols:
raise SymbolsError("Redundant symbols were given")
@@ -252,7 +291,7 @@ def _try_heuristics(f):
result += roots_quadratic(f)
elif n == 3 and flags.get('cubics', True):
result += roots_cubic(f)
- elif n == 4 and flags.get('quartics', False):
+ elif n == 4 and flags.get('quartics', True):
result += roots_quartic(f)
return result
diff --git a/sympy/polys/tests/test_polynomial.py b/sympy/polys/tests/
test_polynomial.py
index 7989bb6..c00216e 100644
--- a/sympy/polys/tests/test_polynomial.py
+++ b/sympy/polys/tests/test_polynomial.py
@@ -1004,27 +1004,52 @@ def test_roots_quadratic():
[-1 + I*sqrt(2)/2, -1 - I*sqrt(2)/2]
def test_roots_cubic():
- assert roots_cubic(Poly(2*x**3, x)) == [0, 0, 0]
- assert roots_cubic(Poly(x**3-3*x**2+3*x-1, x)) == [1, 1, 1]
-
- assert roots_cubic(Poly(x**3+1, x)) == \
- [-1, S.Half - I*sqrt(3)/2, S.Half + I*sqrt(3)/2]
+ #normalized cubic, x**3+a*x**2+b*x+c, is changed to y**3+p*y+q
+ #check c==0
+ assert roots_cubic(Poly(x**3+2*x**2+3*x+0, x)) == [-1 +
I*2**S.Half, 0, -1 - I*2**S.Half]
+ #p==0 and q == 0
+ assert roots_cubic(Poly((x+2)**3,x)) == [-2, -2, -2]
+ #p==0 and q!=0
+ assert set(roots_cubic(Poly(x**3+1, x))) == \
+ set([-S.One, S.Half - I*sqrt(3)/2, S.Half + I*sqrt(3)/2])
+ #p!=0 and q==0
+ assert set(roots_cubic(Poly(x**3 + 3*x**2 + 4*x + 2, x))) == \
+ set([-S.One - I, -S.One, -S.One + I])
+ #p!=0 and q!=0
+ p = Poly(x**3 - 6*x**2 + 13*x - 8, x)
+ ans = roots_cubic(p)
+ assert [[str(round(tmp.evalf().as_real_imag()[j],2)) for tmp in
ans]
+ for j in [0,1]] == [['1.0', '2.5', '2.5'],
+ ['0.0', '-1.32', '1.32']]
def test_roots_quartic():
- assert roots_quartic(Poly(x**4, x)) == [0, 0, 0, 0]
- assert roots_quartic(Poly(x**4 + x**3, x)) in [
- [-1,0,0,0],
- [0,-1,0,0],
- [0,0,-1,0],
- [0,0,0,-1]
- ]
- assert roots_quartic(Poly(x**4 - x**3, x)) in [
- [1,0,0,0],
- [0,1,0,0],
- [0,0,1,0],
- [0,0,0,1]
- ]
- assert roots_quartic(Poly(x**4 + x, x)) == [S.Half + I*sqrt(3)/2,
S.Half - I*sqrt(3)/2, 0, -1]
+ assert roots_quartic(Poly(x**4, x)) == [0]*4
+ # x**4+a*x**3+b*x**2+c*x+d=0 reduces to y^4+e*y^2+f*y+g when
+ # replacing x with y-a/4. Check different types of solutions:
+ # check for homogeneous in x (d=0)
+ set(roots_quartic(Poly(x**4 - 6*x**3 + 11*x**2 - 6*x, x))) == \
+ set([0, 1, 2, 3])
+ # check for homogeneous in y (g=0)
+ set(roots_quartic(Poly(x**4 +4*x**3 + 7*x**2 + 8*x + 4, x))) == \
+ set([-1, -2, -S.Half - I*7**S.Half/2, -1/2 + I*7**S.Half/
2])
+ # no real roots and f=0
+ assert roots_quartic(Poly(10 + 14*x + 11*x**2 + 4*x**3 + x**4,
x)) == \
+ [-1 - I, -1 - 2*I, -1 + I, -1 + 2*I]
+ # no real roots, but f!=0
+ p = Poly(10 + 14*x + 11*x**2 + 4*x**3 + x**4, x)
+ assert set(roots_quartic(p)) == \
+ set([-1 - 2*I, -1 + 2*I, -1 - I, -1 + I])
+ # 2 real roots
+ p = Poly(x**4 - 9*x**3 + 30*x**2 - 42*x + 20, x)
+ ans = roots_quartic(p)
+ assert [[round(tmp.evalf().as_real_imag()[j],2) for tmp in ans]
for j in [0,1]] == \
+ [[1.0, 2.0, 3.0, 3.0],
+ [0.0, 0.0, -1.0, 1.0]]
+ # 4 real roots
+ p = Poly(x**4 + -11*x**3 + 41*x**2 - 61*x + 30, x)
+ ans = roots_quartic(p)
+ assert [[round(tmp.evalf().as_real_imag()[j],2) for tmp in ans]
for j in [0,1]] == \
+ [[1.0, 2.0, 3.0, 5.0], [0.0, -0.0, 0.0, -0.0]]
def test_roots_binomial():
assert roots_binomial(Poly(5*x, x)) == [0]
@@ -1113,8 +1138,7 @@ def test_roots():
r1_2, r1_3, r1_9, r4_9, r19_27 = [ Rational(*r) \
for r in ((1,2), (1,3), (1,9), (4,9), (19,27)) ]
- assert roots(x**3+x**2-x+1, x, cubics=True) in [
- {
+ ans = {
-r1_3 - (r19_27 + r1_9*3**r1_2*11**r1_2)**r1_3 - \
r4_9*(r19_27 + r1_9*3**r1_2*11**r1_2)**(-r1_3): 1,
@@ -1124,41 +1148,33 @@ def test_roots():
-r1_3 + r1_2*(r19_27 + r1_9*3**r1_2*11**r1_2)**r1_3 + \
r4_9/(r1_2 - r1_2*I*3**r1_2)*(r19_27 + r1_9*3**r1_2*11**r1_2)
**(-r1_3) - \
- r1_2*I*3**r1_2*(r19_27 + r1_9*3**r1_2*11**r1_2)**r1_3: 1,
- },
- {
- -r1_3 - (r19_27 + r1_9*3**r1_2*11**r1_2)**r1_3 - \
- r4_9*(r19_27 + r1_9*3**r1_2*11**r1_2)**(-r1_3): 1,
-
- -r1_3 + r1_2*(r19_27 + r1_9*3**r1_2*11**r1_2)**r1_3 - \
- r4_9/(-r1_2 - r1_2*I*3**r1_2)*(r19_27 + r1_9*3**r1_2*11**r1_2)
**(-r1_3) + \
- r1_2*I*3**r1_2*(r19_27 + r1_9*3**r1_2*11**r1_2)**r1_3: 1,
+ r1_2*I*3**r1_2*(r19_27 + r1_9*3**r1_2*11**r1_2)**r1_3: 1
+ }
- -r1_3 + r1_2*(r19_27 + r1_9*3**r1_2*11**r1_2)**r1_3 + \
- r4_9/(r1_2 - r1_2*I*3**r1_2)*(r19_27 + r1_9*3**r1_2*11**r1_2)
**(-r1_3) - \
- r1_2*I*3**r1_2*(r19_27 + r1_9*3**r1_2*11**r1_2)**r1_3: 1,
- },
- ]
+ assert set([t.evalf() for t in roots(x**3+x**2-x+1, x,
cubics=True)]) == \
+ set([t.evalf() for t in ans])
f = (x**2+2*x+3).subs(x, 2*x**2 + 3*x).subs(x, 5*x-4)
r1_2, r13_20, r1_100 = [ Rational(*r) \
for r in ((1,2), (13,20), (1,100)) ]
- assert roots(f, x) == {
+ ans = {
r13_20 + r1_100*(25 - 200*I*2**r1_2)**r1_2: 1,
r13_20 - r1_100*(25 - 200*I*2**r1_2)**r1_2: 1,
r13_20 + r1_100*(25 + 200*I*2**r1_2)**r1_2: 1,
r13_20 - r1_100*(25 + 200*I*2**r1_2)**r1_2: 1,
}
+ assert set([t.evalf() for t in roots(f, x)]) == \
+ set([t.evalf() for t in ans])
p = Poly(z**3 + (-2 - y)*z**2 + (1 + 2*y - 2*x**2)*z - y +
2*x**2, z)
- assert roots(p) == {
+ assert set((t.evalf() for t in roots(p))) == set((t.evalf() for t
in {
S.One: 1,
S.Half + S.Half*y + S.Half*(1 - 2*y + y**2 + 8*x**2)**S.Half:
1,
S.Half + S.Half*y - S.Half*(1 - 2*y + y**2 + 8*x**2)**S.Half:
1,
- }
+ }))
assert roots(a*b*c*x**3 + 2*x**2 + 4*x + 8, x, cubics=False) ==
{}
assert roots(a*b*c*x**3 + 2*x**2 + 4*x + 8, x, cubics=True) != {}
diff --git a/sympy/simplify/tests/test_simplify.py b/sympy/simplify/
tests/test_simplify.py
index a8733f3..cca891e 100644
--- a/sympy/simplify/tests/test_simplify.py
+++ b/sympy/simplify/tests/test_simplify.py
@@ -90,6 +90,8 @@ def test_factorial_simplify():
def test_simplify():
x,y,z,k,n,m,w,f,s,A = symbols('xyzknmwfsA')
+ assert all(simplify(tmp)==tmp for tmp in [I,E,oo,x,-x,-oo,-E,-I])
+
e = 1/x + 1/y
assert e != (x+y)/(x*y)
assert simplify(e) == (x+y)/(x*y)
--
1.6.3.2.1299.gee46c
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"sympy-patches" group.
To post to this group, send email to sympy-patches@googlegroups.com
To unsubscribe from this group, send email to
sympy-patches+unsubscr...@googlegroups.com
For more options, visit this group at
http://groups.google.com/group/sympy-patches?hl=en
-~----------~----~----~----~------~----~------~--~---
