New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Sorting of polynomial roots #14293

Open
skirpichev opened this Issue Feb 20, 2018 · 57 comments

Comments

Projects
None yet
4 participants
@skirpichev
Contributor

skirpichev commented Feb 20, 2018

In fact, it's about their identification.

Low-level code in rootisolation module produces complex roots, sorted by coordinates of south-western corner of the isolation rectangle (by ax with ties broken by ay). Clearly, this sorting is not stable (in sense that refining process doesn't preserve one) => no way to index roots. Yet, this sorting has some nice properties for polynomials over Q, which are preserved, namely - complex roots coming as conjugated pairs, making conjugate() a trivial task.

Incredibly cumbersome code in RootOf, introduced by (unreviewed commit) 73fafba - makes situation only worse. Apparently, it doesn't solve stability sorting issue (Chris, the problem is not just with pure imaginary roots!), introduces bugs (#14291) and breaks conjugated pairs, so they are not adjacent anymore, e.g.:

In [1]: Poly((x**2 - 2)*(x**2 + 2)*(x**2 + 3)).all_roots()
Out[1]: [-√2, √2, -√3⋅ⅈ, -√2⋅ⅈ, √2⋅ⅈ, √3⋅ⅈ]

Is there some algorithm that sort roots in a some stable way? If so, it should be used to index roots. If no, individual roots should be specified differently, e.g. by isolation rectangle. Not by index. Perhaps, a third option would be ordering by ax and handle rest symbolically.

@smichr

This comment has been minimized.

Member

smichr commented Feb 21, 2018

Perhaps sorting by (x, abs(y), sign(y)) would be preferable. ordered can easily and efficiently use such a key.

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Feb 21, 2018

@smichr, please explain what do you mean by x, y. E.g., if south-western corner of the isolation rectangle (real and imaginary parts, respectively), then this will be not better that sorting used currently by rootisolation module.

Let me recall. The issue is: sorting should allow you to index roots. refine() of isolation rectangles shouldn't break this order.

@smichr

This comment has been minimized.

Member

smichr commented Feb 21, 2018

sorting should allow you to index roots. refine() of isolation rectangles shouldn't break this order.

I'm not sure I understand what it is the most important issue to you. Here are my assumptions:

  1. some polynomial p(x) has roots.
  2. indexing is arbitrary and is not a mathematical property of the root
  3. it might make sense to make indexing canonical
  4. you don't know what a roots index will be until you define the rules
  5. once you have rules you won't know what the index is until you determine the value of the root
  6. values need to be computed with enough precision so they can be unambiguously distinguished from others; this is the process of refinement
  7. a refined root can be assigned an index according to the rules; an unrefined root cannot be assigned an index

example

>>> sorted(set([randint(-3,3) + I*randint(-3,3) for i in range(10)]),
... key=lambda x: (abs(im(x))!=0, abs(re(x))!=0, re(x), abs(im(x)), sign(im(x))))
[−1,3,−3i,3i,−3−2i,−2−i,−1+i,−1−2i,1+2i,2−3i]
@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Feb 21, 2018

it might make sense to make indexing canonical

We want this, as currently it's the only way to specify individual roots via RootOf.

The question: is there an algorithm, that allow us to "unambiguously distinguish" (after some finite number of refinement process steps) roots to enumerate them (read "assign index"), in order that will be unchanged by further refinement of roots.

a refined root can be assigned an index according to the rules; an unrefined root cannot be assigned an index

That might be an option, i.e. partially assign indexes and handle rest of roots symbolically (i.e. no ordering, no evalf, etc). But that more or less kills the purpose of index... - then intervals must be used by RootOf instead as primary syntax to distinguish roots.

@smichr

This comment has been minimized.

Member

smichr commented Feb 22, 2018

is there an algorithm, that allow us to "unambiguously distinguish"

Please correct me if I am wrong, but I believe the current method is such an algorithm. Roots are refined until they are distinct from one another and then are assigned indices. Are you wanting to assign an index before doing any refinement?

breaks conjugated pairs, so they are not adjacent anymore

Yes, the reals are in order from negative to positive, the imaginaries are in order from negative to positive magnitude and the rest are in pairs.

>>> print(Poly((x+I)*(x-I)*(x+1)*(x-3)*(x-2)*(2*x-3*I)*(2*x+3*I)*(x-(3+4*I))*(x-(3+4*I).conjugate())*(x-(1+4*I))*(x-(1+4*I).conjugate())).all_roots())
[-1, 2, 3, -3*I/2, -I, I, 3*I/2, 1 - 4*I, 1 + 4*I, 3 - 4*I, 3 + 4*I]

I take it that you would like to have the property be that index i and i + 1 are the conjugate pairs for i in range(nreal, len(allroots)?

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Feb 22, 2018

Please correct me if I am wrong, but I believe the current method is such an algorithm.

Oh. The word "algorithm" means (in particular), that it will end computation after several (finite number!) steps. It can't stuck like #14291.

Roots are refined until they are distinct from one another

Could you formalize here criteria for "being distinct", e.g. in terms of isolating rectangles? Apparently, if there are some algorithm in your code - it was wrongly implemented.

Are you wanting to assign an index before doing any refinement?

I want to assign indexes after some well defined procedure of refinement. Which will have a defined rules of termination and will terminate in a finite number of steps.

I take it that you would like to have the property be that index i and i + 1 are the conjugate pairs

The whole issue is not about this, but yes. Such ordering (which used by low-level rootisolation module) has some benefits.

@smichr

This comment has been minimized.

Member

smichr commented Feb 23, 2018

Ok, I understand the issue. No solution, yet.

@smichr

This comment has been minimized.

Member

smichr commented Feb 23, 2018

Given irreducible polynomial P(x) with real coefficients we know that all non-real roots come in pairs.

If two pairs share the same real or imaginary coordinate then P(x) is not irreducible:

>>> c, d, e = var('c:e', real=True)
>>> eq
(-a + x)*(-b + x)*(x - conjugate(a))*(x - conjugate(b))
>>> factor(eq.subs(a, c+e*I).subs(b,c+d*I).expand())
(c**2 - 2*c*x + d**2 + x**2)*(c**2 - 2*c*x + e**2 + x**2)
>>> factor(eq.subs(a, d+c*I).subs(b,e+c*I).expand())
(c**2 + d**2 - 2*d*x + x**2)*(c**2 + e**2 - 2*e*x + x**2)

Since none of the pairs will share the same real or imaginary value the non-real complex roots of an irreducible polynomial can be refined numerically into non-intersecting rectangles.

So at least a first step is to not process the roots of more than one P(x) at a time. Sort roots by the expression that they come from and then sub-sort them by the refined values, arbitrarily using the imaginary values, putting negative before positive after identifying the conjugate pairs (which have overlapping x ranges for their bounding rectangle).

@asmeurer

This comment has been minimized.

Member

asmeurer commented Feb 24, 2018

@smichr I'm not sure if your proof is correct. A factorization of a polynomial only shows that it is irreducible (over Q) if the coefficients are rational. But there's no reason that c**2 + d**2 and so on should be rational. With that being said, I don't know if the statement is true or not (except for the trivial case of a double real root, which is easy to show by the square free decomposition).

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Feb 24, 2018

Given irreducible polynomial P(x) with real coefficients we know that all non-real roots come in pairs.

That's true for any polynomial with real coefficients.

If two pairs share the same real or imaginary coordinate then P(x) is not irreducible:

As I said, your statement above is valid irrespectively of irreducibility if we consider our polynomial over the field of coefficients big enough, namely, over C. Therefore, this statement already wrong.

But it's easy to construct a counterexample:

In [7]: x**4 + 10*x**2 + 1
Out[7]: 
 4       2    
x  + 10⋅x  + 1

In [8]: factor(_)
Out[8]: 
 4       2    
x  + 10⋅x  + 1

In [9]: factor(_7, extension=sqrt(2) + sqrt(3))
Out[9]: 
⎛ 2       ___    ⎞ ⎛ 2       ___    ⎞
⎝x  - 2⋅╲╱ 6  + 5⎠⋅⎝x  + 2⋅╲╱ 6  + 5⎠

In [10]: roots(_7)
Out[10]: 
⎧      _______________          _______________           _____________          _____________   ⎫
⎨     ╱       ___              ╱       ___               ╱     ___              ╱     ___        ⎬
⎩-ⅈ⋅╲╱  - 2⋅╲╱ 6  + 5 : 1, ⅈ⋅╲╱  - 2⋅╲╱ 6  + 5 : 1, -ⅈ⋅╲╱  2⋅╲╱ 6  + 5 : 1, ⅈ⋅╲╱  2⋅╲╱ 6  + 5 : 1

There are two conjugated pairs, that share same real part.

@smichr

This comment has been minimized.

Member

smichr commented Feb 26, 2018

Can you construct an example in which the real part is not 0?

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Feb 27, 2018

Sure. Shift this polynomial by real number you want, isn't this clear?

@smichr

This comment has been minimized.

Member

smichr commented Feb 27, 2018

isn't this clear?

Yes, and I see the problem. So is there any way to know that a polynomial for which factor(p(x)).is_Add is true has real and imaginary parts that are the same (even though you might not know what they are)?

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Feb 27, 2018

I hope, we shouldn't care about equal imaginary parts.

@smichr

This comment has been minimized.

Member

smichr commented Feb 27, 2018

shouldn't care about equal imaginary parts.

If that is the case (and I'm not sure why) then perhaps the answer is clear in terms of identifying roots: refine until either the real or imaginary parts are distinct for each expression.

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Feb 27, 2018

If that is the case (and I'm not sure why)

It is.

refine until either the real or imaginary parts are distinct for each expression.

I told in the issue description (I hope this time it will be clear), that such ordering (which is used by rootisolation module) is not stable. If you continue refinement process, this order may change. So, you can't index roots this way.

@smichr

This comment has been minimized.

Member

smichr commented Feb 28, 2018

,

issue description (I hope this time it will be clear), that such ordering (which is used by rootisolation module) is not stable

PR #14348 attempts to provide some level of stability in that a root and its conjugate are always next to each other with the one having the negative imaginary part having the lower index. An example of any polynomial whose indices are not stable wrt refinement would be appreciated feedback...I should know by tomorrow if failing tests provide such feedback.

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Feb 28, 2018

PR #14348 attempts to provide some level of stability in that a root and its conjugate are always next to each other

For this you could just revert your patch 73fafba. The rootisolation module already does such sorting.

An example of any polynomial whose indices are not stable wrt refinement would be appreciated feedback...

Could you specify how it's different from the sorting in rootisolation.py and why it's different?

skirpichev added a commit to skirpichev/diofant that referenced this issue Feb 28, 2018

This reverts 73fafba
See sympy/sympy#14293

There is no stable sorting yet, further refinement **may**
change order of roots.  Setting ``eps`` only hide
problem a little.
@smichr

This comment has been minimized.

Member

smichr commented Feb 28, 2018

The structure of the roots(p(x)) is now [sorted real roots, complex pairs of f1(x), complex pairs of f2(x), ...] where p(x) = f1(x)*f2(x)....

Reverting patch 73fafba will not produce this sort of ordering; the tests added in that patch were the motivation at that time. If you revert the patch, try them or any of the ones raised in this issue to see the problem again.

>>> a=Poly((x**2 - 2)*(x**2 + 2)*(x**2 + 3))
>>> b=Poly(x**4 + 10*x**2 + 1)
>>> c=Poly(((x - 1)**2 + 1)*((x - 1)**2 + 2)*(x - 1))
>>> [i.n(3, chop=True) for i in (a*b*c).all_roots()]
[-1.41, 1.00, 1.41, 
-1.41*I, 1.41*I, 
-1.73*I, 1.73*I, 
1.0 - 1.0*I, 1.0 + 1.0*I, 
1.0 - 1.41*I, 1.0 + 1.41*I, 
-0.318*I, 0.318*I, 
-3.15*I, 3.15*I]
>>>
@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Feb 28, 2018

Reverting patch 73fafba will not produce this sort of ordering

I see. But you could just slightly adapt sorting key: (default_sort_key(r[1]), r[0].ax, r[0].ay) or something like that.

the tests added in that patch were the motivation at that time.

I'm a bit unsure, that this was a motivation. Anyway, another "solution" - better refine initial rectangles, using eps. It seems, with eps=1/10 - current tests pass.

@smichr

This comment has been minimized.

Member

smichr commented Feb 28, 2018

better refine initial rectangles, using eps

I don't see a reason to try refine non-real or non-imaginary roots whose sorting will inevitably be arbitrary. Once they have been refined enough to be uniquely identified, your work is done

I suspect, also, that you would easily come up with a case where eps=1/10 was not sufficient to distinguish between two roots. The is_disjoint (sp?) tests is the only way to know that the intervals are refined enough and the algorithm currently is set up to refine until that criteria is affirmed for each pair of roots.

I don't see that any sorting key that does not separate the roots according to underlying expression from which they came will succeed at yielding a stable sort. If two ranges overlapped in reals but were distinct in imaginary ranges they would be disjoint -- but the indices of the roots might flip with further refinement. I don't believe this flip can happen if the roots are sorted by subexpression/factor: roots of a given factor will either have the same real ranges and easily be identified as conjugate pairs or else the real ranges will be easily made disjoint. There will be no flipping of indices for that expression at greater refinement. And I think this is the property you identified as being desirable in the issue. In addition, changes made in the PR I propose intend to guarantee that the conjugate pairs are always next to each other -- the _eval_conjugate method has been added to RootOf to take advantage of that property so, for example, conjugate(RootOf(e, i)) will always be RootOf(e, i + 1) if the imaginary part of the conjugate argument is negative else RootOf(e, i - 1).

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Feb 28, 2018

I don't see a reason to try refine non-real or non-imaginary roots whose sorting will inevitably be arbitrary.

The reason is not break tests too much.

I suspect, also, that you would easily come up with a case where eps=1/10 was not sufficient to distinguish between two roots.

Sure, as I said, this variant is not a real solution.

I don't see that any sorting key that does not separate the roots according to underlying expression from which they came will succeed at yielding a stable sort.

Variant with a different sort key solves same problem as your PR (i.e., "the structure of the roots(p(x)) is now [sorted real roots, complex pairs of f1(x), complex pairs of f2(x), ...]"), but it costs just a few lines.

I don't believe this flip can happen if the roots are sorted by subexpression/factor

Why not? It can be. "Flips" will occur not only between roots from same factors, but between different factors too. The reason is that different factors may share close roots, no problem.

changes made in the PR I propose intend to guarantee that the conjugate pairs are always next to each other

This is the only property, that can validate your changes. But it also could be preserved just by different sort key (sort by factor, then by ax, then by ay).

@smichr

This comment has been minimized.

Member

smichr commented Feb 28, 2018

I don't think that the variant sort key, without a refinement process, is sufficient to assure that further refinement will not cause a flip.

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Mar 1, 2018

I don't think that the variant sort key, without a refinement process, is sufficient to assure that further refinement will not cause a flip.

Neither solution, which I know, allow you avoid such "flips".

@smichr

This comment has been minimized.

Member

smichr commented Mar 1, 2018

If you mean by flip that one cannot assign an index and have it remain invariant for a desired sorting of the roots (like key=lambda r: (re(r), im(r))) then I agree. But I don't think the task of assigning a unique ordering is difficult. Is there a problem with the following? It is very close to what SymPy is doing without the offending commit.

  1. P(x) is irreducible; none of its roots are repeated
  2. roots of P are real or complex
  3. the real roots can be sorted on the number line
  4. the non-real roots have an arbitrary ordering
  5. the non-real roots occur in pairs
  6. let n be the degree of P and h = n/2
  7. there is a finite extend (rectangle) of the lower half of the complex plane that will contain h roots
  8. use a binary division of the rectangle to isolate the roots
  9. stop dividing once the h roots have been isolated into non-intersecting rectangles
  10. label the furthest rectangle to the left as 0
  11. label all rectangles that overlap with the horizontal bounds of rectangle 0 as 1..n in increasing distance from the x-axis
  12. continue with the n+1st rectangle and all that overlap with it, etc...
  13. when you are done you will have rectangles numbered 0 through h-1
  14. those numbers are arbitrary but they are assigned deterministically
  15. those numbers j = 0, 1, ..., h - 1 correspond to non-real roots with negative imaginary parts
  16. the conjugate roots in the upper half plane have bounding rectangles given by the reflection of the lower rectangles across the x axis
  17. number the conjugate pairs (j_lower, j_upper) as (2j, 2j + 1)

The numbers that are obtained should be determined by the bounds chosen at the outset and the method used to subdivide that space. Together, those constitute an "algorithm" to uniquely number the roots of P. It guarantees nothing more than that conjugates will be next to each other.

skirpichev added a commit to skirpichev/diofant that referenced this issue Mar 1, 2018

This reverts 73fafba
See sympy/sympy#14293

There is no stable sorting yet, further refinement **may**
change order of roots.  Setting ``eps`` only hide
problem a little.
@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Mar 1, 2018

Yes, it seems in this way roots of irreducible, square-free polynomial can be numbered. More or less on same way it may work for polynomials with complex coefficients (e.g. gaussian rationals).

But this numbering depends on initial bounds chosen. Initial rectangle should be generated by some specific scheme, not just include all roots.

For reducible polynomials, probably, you can extend this "method", sorting roots by factor and then by root index "inside" of this factor.

But I think, that there is a way to sort roots in a mathematically meaningful way. After all, in principle we can solve zero-decision problem for algebraic numbers. So, it should be decidable if real parts of roots are same or not.

skirpichev added a commit to skirpichev/diofant that referenced this issue Mar 1, 2018

This reverts 73fafba
See sympy/sympy#14293

There is no stable sorting yet, further refinement **may**
change order of roots.  Setting ``eps`` only hide
problem a little.
@smichr

This comment has been minimized.

Member

smichr commented Mar 27, 2018

you decide somehow if they share same real part

No. I never mix roots of different factors and I separate real and imaginary from those that are neither.

  • reals are solved without problem
  • imaginary roots are solved as real roots by working with poly(I*x) instead of poly(x)
  • complex roots can be refined with assurance that the root bounds will resolve into distinct quadrilaterals, e.g.
>>> [i.n(2) for i in Poly((x**4 + 10*x**2 + 1).subs(x,x-1)).all_roots()]
[1.0 - 0.32*I, 1.0 + 0.32*I, 1.0 - 3.1*I, 1.0 + 3.1*I]
@smichr

This comment has been minimized.

Member

smichr commented Mar 27, 2018

Also, roots are orderly though not sorted as noted.

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Mar 27, 2018

Ok, the issue is not solved yet.

@jksuom

This comment has been minimized.

Member

jksuom commented Mar 27, 2018

I'm not sure how multiplicity of real roots of equation for u is connected with the original equation.

It is a classical problem in intersection theory that the number of solutions cannot be computed directly from the equations obtained by elimination. (The solution was "moving the axes to general position" so that for each solution of the eliminated equations there is only one solution of the original equations.)

But the number of complex solutions having the same real part x can be computed otherwise. Take a (rational) interval [a, b] containing x but no other real part of any solution and compute the number of roots in a long quadrilateral [a, b] x [-M, M] where M is an upper bound for the absolute value of all roots.

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Mar 27, 2018

Ok, so we can just ignore multiplicity of real roots equation for u.

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Mar 27, 2018

So, finally, following should solve the sorting problem:

  1. Compute isolation intervals for real parts of every factor.
  2. Refine them to make disjoint (exactly as we do for real roots now).
  3. Then refine complex intervals to be both (a) pairwise disjoint (as now) and (b) real part of every such interval is inside of exactly one interval, computed on steps 1-2.

@jksuom, did I miss something?

However, there may be some performance penalties if we will run these steps instead of current sorting: computing isolation intervals for real parts in the way you suggest will produce polynomials of high degree (e.g. in this example r has degree 49).

@jksuom

This comment has been minimized.

Member

jksuom commented Mar 28, 2018

That looks like what I had imagined. Performance will probably suffer but isolating the real roots should not be too expensive. (I would expect Sturmian methods to be fairly efficient. I would also consider looking into subresultants_qq_zz.)

@smichr

This comment has been minimized.

Member

smichr commented Mar 28, 2018

will produce polynomials of high degree

Did you take a look at the mathstackexchange discussion? Would the following be advantageous in not having to solve such high order polynomials (because you would only need to solve the gcd)?

def samere(p, r):
  def res(p):
    a, b = var('a b', real=True)
    p1a, p1b = p.subs(x, a + I*b).as_real_imag()
    return Poly(p1a, b).resultant(Poly(p1b, b))
  g = Poly(res(p)).gcd(Poly(res(r)))
  if g.degree() > 0:
    return set(g.real_roots())


>>> p
x**4 - 8*x**3 + 26*x**2 - 40*x + 13
>>> r
x**4 - 8*x**3 + 30*x**2 - 56*x + 29
>>> samere(p, r)
{2}
>>> [i.n(2) for i in Poly(p).all_roots()]
[0.43, 3.6, 2.0 - 2.1*I, 2.0 + 2.1*I]
>>> [i.n(2) for i in Poly(r).all_roots()]
[0.79, 3.2, 2.0 - 2.7*I, 2.0 + 2.7*I]

I would imagine the sorting procedure to be:

  1. have roots that have been made disjoint for each factor (as in my current PR)
  2. sort real and imaginary
  3. for any pair of non-real, non-imag roots (from different factors) that are not disjoint, see if they share a common real root that falls within their real bounds and
    a) if so, they are a tie and are refined until their imaginary parts are different while setting their real part to be the intersection of the current bound, e.g. (1, 3) & (2,4) becomes (2, 3); or
    b) if not, refine them each until they are distinct in real parts
  4. sort by (ivl.ax, ivl.ay)
  5. merge imaginary into non-real/non-imaginary to give non-real roots
  6. return reals + non-real
@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Mar 28, 2018

Did you take a look at the mathstackexchange discussion?

Sure.

Would the following be advantageous in not having to solve such high order polynomials

Of course, no. gcd couldn't reduce degree of resultant for same factor. E.g. x**7 + I*x**4 + 2*x - I from my example is irreducible.

have roots that have been made disjoint for each factor (as in my current PR)

They are already coming as disjoint from rootisolation routines, if being "disjoint" means have pairwise non-intersecting rectangles.

@smichr

This comment has been minimized.

Member

smichr commented Mar 28, 2018

already coming as disjoint from rootisolation routines, if being "disjoint" means have pairwise non-intersecting rectangles.

I don't believe this is true in the proposed PR: they are only pairwise disjoint for a given factor. Yes, all roots are disjoint.

@smichr

This comment has been minimized.

Member

smichr commented Mar 30, 2018

So here is a proof of concept showing how I would sort roots in the same order as nsort (reals then by (re, im) for complex roots). I tested it against several 100 polynomials without failures (other than those due to rounding errors). It takes about 8X longer to sort than it does to compute the roots.

See here.

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Mar 30, 2018

It seems, you still think that for same poly "normal refinement will work"... This makes me sad.

I tested it against several 100 polynomials without failures (other than those due to rounding errors).

Oh. Are you sure that even all lines are covered?

BTW, you must test mathematical ordering, not nroots. And there should be no errors.

@smichr

This comment has been minimized.

Member

smichr commented Mar 30, 2018

Thanks for giving this a close look. And yes, it should have been _nsort not nroots though it looks like the two are giving the same result. Also, I agree that there should be no (algorithm related) errors -- which I was able to confirm that there are. (The errors I saw before were rounding errors.)

still think that for same poly "normal refinement will work"

Although I may have let the assumption slip through it is not my intention to have done so. I understand that it can be expected to succeed when two intervals do not share a common real root. Any attempt to do otherwise is simply a shortcoming of my rather limited mathematical insight. "I'm a doctor, not a [mathematician]!" :-)

I'll take a look again in a couple weeks. Current state is here.

May any who read find renewed hope in what was won, one third day, so many years ago.

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Mar 30, 2018

it should have been _nsort not nroots though it looks like the two are giving the same result.

I did mention that because it's not clear whether "failures due to rounding" are connected with nroots and not with mistakes in your "algorithm".

Current state is here.

I can't run this.

In [1]: from a import rsort

In [2]: rsort(Poly(-3*x**5 - 3*x**4 - 3*x**3 + x**2 + 2*x + 1).all_roots())
~/src/sympy/sympy/core/numbers.py in __new__(cls, p, q, gcd)
   1510 
   1511             if not isinstance(p, SYMPY_INTS + (Rational,)):
-> 1512                 raise TypeError('invalid input: %s' % p)
   1513             q = q or S.One
   1514             gcd = 1

TypeError: invalid input: 0

BTW, I don't see if your code does now something meaningful, even in case p1a!=p1b (different factors). Ok, you take first real part you have found, that match both projections of complex root rectangles to the real axis. Fine. But what if there are another such real part?

@smichr

This comment has been minimized.

Member

smichr commented Mar 31, 2018

If there is more than 1 root in the overlapping ranges, just refine until there is 1...then set the overlapping root's xranges to the intersection of the overlapping x ranges.

@smichr

This comment has been minimized.

Member

smichr commented Mar 31, 2018

I have not implemented that root xrange update in the posted gist. I was running that code in the pr branch of mine named cr...that allows Rational to handle the mptypes without using str.

skirpichev added a commit to skirpichev/diofant that referenced this issue Mar 31, 2018

skirpichev added a commit to skirpichev/diofant that referenced this issue Apr 5, 2018

@smichr

This comment has been minimized.

Member

smichr commented Apr 11, 2018

Here's a test suite. (I tried to download diofant to try it, but ran into too many difficulties.)

def test_rsort():
    ps = (59*x**5 - 35*x**4 + 56*x**3 + 39*x**2 - 23*x + 27,
        64*x**5 - 17*x**4 + 72*x**3 + 27*x**2 - 3*x + 10,
        74*x**5 - 33*x**4 + 99*x**3 + 43*x**2 - 2*x + 42,
        -3*x**5 - 3*x**4 - 3*x**3 + x**2 + 2*x + 1,
        61*x**5 - 20*x**4 + 34*x**3 - 100*x**2 - 40*x - 18,
        50*x**5 + 19*x**4 + 63*x**3 - 63*x**2 - 41*x - 40,
        x**7 - x**6 + 10*x**5 - 7*x**4 + 22*x**3 - 10*x**2 - 13*x - 2,
        x**8 - 16*x**7 + 120*x**6 - 544*x**5 + 1590*x**4 - 2992*x**3 +
        3384*x**2 - 1888*x + 377)
    for p in ps:
        R = Poly(p).all_roots()
        a = [str(i.n(4)) for i in (R)]
        b = [str(i) for i in sorted([i.n(4) for i in R],
            key=lambda x: (
                not x.is_real,
                x.as_real_imag()[0],
                abs(x.as_real_imag()[1]),
                sign(x.as_real_imag()[1])))]
        assert a == b
@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Apr 11, 2018

Here's a test suite.

Hmm, what exactly it tests?

In [30]: roots4 = Poly(x**4 + 10*x**2 + 1).all_roots()

In [31]: roots42 = sorted(roots4, key=lambda x: (
    ...:                                 not x.is_real,
    ...:                                 x.as_real_imag()[0],
    ...:                                 abs(x.as_real_imag()[1]),
    ...:                                 sign(x.as_real_imag()[1])))
    ...:                                 

In [32]: [_.n(2) for _ in roots4]
Out[32]: [1.6e-14 - 3.1⋅ⅈ, 2.6e-26 - 0.32⋅ⅈ, 2.6e-26 + 0.32⋅ⅈ, 1.6e-14 + 3.1⋅ⅈ]

In [33]: [_.n(2) for _ in roots42]
Out[33]: [1.6e-14 + 3.1⋅ⅈ, 2.6e-26 + 0.32⋅ⅈ, 2.6e-26 - 0.32⋅ⅈ, 1.6e-14 - 3.1⋅ⅈ]

(this was in diofant, but I don't think it does matter for the last line)

I tried to download diofant to try it, but ran into too many difficulties.

Works for me (are you trying python 2?), but currently diofant uses different sorting key, so some tests will "fail".

@smichr

This comment has been minimized.

Member

smichr commented Apr 13, 2018

It checks that the returned roots (in whatever branch of yours or mine that sorts the RootOf instance) are sorted the same way that the numerical roots would be sorted with the given key.

>>> roots4 = [i.n(4) for i in Poly(x**4 + 10*x**2 + 1).all_roots()]
>>> roots42 = sorted(roots4, key=lambda x: (
...                                      not x.is_real,
...                                      x.as_real_imag()[0],
...                                      abs(x.as_real_imag()[1]),
...                                      sign(x.as_real_imag()[1])))
>>> roots4
[-0.3178*I, 0.3178*I, -3.146*I, 3.146*I]
>>> roots42
[-0.3178*I, 0.3178*I, -3.146*I, 3.146*I]

diofant appears to need a cache-related library which I didn't have and something about how I tried to install it got me stuck with a git error that required a reboot to escape from.

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented Apr 13, 2018

are sorted the same way that the numerical roots would be sorted with the given key.

It appears to be a mathematical sorting. Do you believe you have solved the issue?

roots4
[-0.3178I, 0.3178I, -3.146I, 3.146I]

Is this on your branch? I got TypeError on cr branch and nothing on the sympy master. Please also check shifted polynomial, it seems you have added some special handling for imaginary roots.

something about how I tried to install it got me stuck with a git error

Have you tried to follow installation instructions?

@smichr

This comment has been minimized.

Member

smichr commented Apr 14, 2018

Do you believe you have solved the issue?

No. The cr branch does not involve sorting. I was working on a sorting routine, however, and that's where this test came from. If you are getting a TypeError, perhaps you are not using the most up-to-date version of the cr branch.

check shifted
Like this?

>>> [i.n(4) for i in Poly((x**4 + 10*x**2 + 1).subs(x,x-1)).all_roots()]
[1.0 - 0.3178*I, 1.0 + 0.3178*I, 1.0 - 3.146*I, 1.0 + 3.146*I]

Have you tried to follow installation instructions?

I was able (now) to follow them. On Windows (for developer version) I had to do the first step in my git program and the 2nd step worked in a command window but not from within the git program shell:

MINGW32 ~/diofant (master)
$ pip install -e .[develop,docs]
bash: pip: command not found

MINGW32 ~/diofant (master)
$ python pip install -e .[develop,docs]
C:\Python36-32\python.exe: can't open file 'pip': [Errno 2] No such file or directory

I was then able to run everything as expected.

>>> [i.n(4,chop=True) for i in Poly(x**4 + 10*x**2 + 1).all_roots()]
[-3.146*I, -0.3178*I, 0.3178*I, 3.146*I]

skirpichev added a commit to skirpichev/diofant that referenced this issue May 6, 2018

skirpichev added a commit to skirpichev/diofant that referenced this issue May 15, 2018

skirpichev added a commit to skirpichev/diofant that referenced this issue May 18, 2018

skirpichev added a commit to skirpichev/diofant that referenced this issue May 20, 2018

@skirpichev

This comment has been minimized.

Contributor

skirpichev commented May 24, 2018

FYI, there is an attempt to fix this issue (and also #14738) in diofant/diofant#633.

skirpichev added a commit to skirpichev/diofant that referenced this issue May 28, 2018

skirpichev added a commit to skirpichev/diofant that referenced this issue Jun 3, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment