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 -~----------~----~----~----~------~----~------~--~---