# Ticket #10803(closed defect: fixed)

Opened 2 years ago

## critical bug in real_roots

Reported by: Owned by: zimmerma jason, jkantor critical sage-5.4.1 numerical real roots, sd40.5 cwitty, was N/A Jeroen Demeyer Paul Zimmermann sage-5.4.1.rc0

```sage: from sage.rings.polynomial.real_roots import real_roots
sage: x = polygen(QQ)
sage: f = 2503841067*x^13 - 15465014877*x^12 + 37514382885*x^11 - 44333754994*x^10 + 24138665092*x^9 - 2059014842*x^8 - 3197810701*x^7 + 803983752*x^6 + 123767204*x^5 - 26596986*x^4 - 2327140*x^3 + 75923*x^2 + 7174*x + 102
sage: len(real_roots(f)), len(real_roots(f,strategy='warp'))
(11, 13)
```

It seems different seeds give different numbers of roots:

```sage: len(real_roots(f,seed=1))
11
sage: len(real_roots(f,seed=3))
13
```

## Change History

### comment:1 follow-up: ↓ 2 Changed 2 years ago by dsm

Note that it's not just that real_roots hasn't fully resolved the roots, and so an interval contains more than one root with the right multiplicity being reported, which would also lower len(real_roots(f)). In the above case two roots are simply missing, as can be confirmed by looking at the reported roots themselves.

### comment:2 in reply to: ↑ 1 Changed 2 years ago by zimmerma

Note that it's not just that real_roots hasn't fully resolved the roots, and so an interval contains more than one root with the right multiplicity being reported, which would also lower len(real_roots(f)). In the above case two roots are simply missing, as can be confirmed by looking at the reported roots themselves.

good point. I should have written:

```sage: sum(x[1] for x in real_roots(f))
11
sage: sum(x[1] for x in real_roots(f,seed=2))
13
```

Paul

### comment:3 Changed 2 years ago by jdemeyer

• Description modified (diff)

### comment:4 Changed 2 years ago by jdemeyer

Let me point out that f.real_roots() does seem to work (and uses a different algorithm).

### comment:5 Changed 2 years ago by mhampton

Here's what might be a clue: somehow calling the roots method on the "ocean" class fixes the problem (I don't yet have any idea why):

```from sage.rings.polynomial.real_roots import *
x = polygen(QQ)
f = 2503841067*x^13 - 15465014877*x^12 + 37514382885*x^11 - 44333754994*x^10 + 24138665092*x^9 - 2059014842*x^8 - 3197810701*x^7 + 803983752*x^6 + 123767204*x^5 - 26596986*x^4 - 2327140*x^3 + 75923*x^2 + 7174*x + 102
```
```ctx = context(False, 1, 32)
b, other = to_bernstein(f, -1, 7)
oc1 = ocean(ctx, bernstein_polynomial_factory_ratlist(b), linear_map(-1, 7))
oc1.find_roots()
len(oc1.roots())
```

outputs 11.

```oc = ocean(ctx, bernstein_polynomial_factory_ratlist(b), linear_map(-1, 7))
_ =  oc.roots()
oc.find_roots()
print len(oc.roots())
```

outputs 13.

### comment:6 Changed 2 years ago by mhampton

Sorry, that was misleading. Its not calling roots(), its the re-use of the context object that seems to be fixing things.

### comment:7 Changed 2 years ago by mhampton

Presumably the re-use of the context object means the seed is not the same, so I guess that's a red herring - my apologies.

### comment:8 follow-up: ↓ 9 Changed 2 years ago by dsm

FYI, here's the smallest case I know:

```sage: f = 535944*x^12 - 3312356*x^11 + 8042902*x^10 - 9521317*x^9 + 5204589*x^8 - 461471*x^7 - 682572*x^6 + 174775*x^5 + 25797*x^4 - 5793*x^3 - 473*x^2 + 17*x + 1
sage: sum(r[1] for r in real_roots(f,seed=1))
10
sage: len(f.real_roots())
12
```

It seems pretty sensitive to the initial bounds: (-1,7) fails, as does (-1,7+1/2^k) for k=13,14,15,17,21,.. etc. Unfortunately debugging is incredibly slow as real_roots.pyx takes forever to compile for me. :-(

### comment:9 in reply to: ↑ 8 Changed 2 years ago by drkirkby

FYI, here's the smallest case I know:

```sage: f = 535944*x^12 - 3312356*x^11 + 8042902*x^10 - 9521317*x^9 + 5204589*x^8 - 461471*x^7 - 682572*x^6 + 174775*x^5 + 25797*x^4 - 5793*x^3 - 473*x^2 + 17*x + 1
sage: sum(r[1] for r in real_roots(f,seed=1))
10
sage: len(f.real_roots())
12
```

FYI,

```Mathematica 7.0 for Sun Solaris x86 (64-bit)

In[1]:= f = 535944*x^12 - 3312356*x^11 + 8042902*x^10 - 9521317*x^9 + 5204589*x^8 - 461471*x^7 - 682572*x^6 + 174775*x^5 + 25797*x^4 - 5793*x^3 - 473*x^2 + 17*x + 1

2         3          4           5           6
Out[1]= 1 + 17 x - 473 x  - 5793 x  + 25797 x  + 174775 x  - 682572 x  -

7            8            9            10            11
>    461471 x  + 5204589 x  - 9521317 x  + 8042902 x   - 3312356 x   +

12
>    535944 x

In[2]:= NSolve[f==0,x]

Out[2]= {{x -> -0.279345}, {x -> -0.169895}, {x -> -0.0836443},

>    {x -> -0.0391184}, {x -> 0.0525278}, {x -> 0.217266}, {x -> 0.608742},

>    {x -> 0.743349}, {x -> 0.955059}, {x -> 0.956605}, {x -> 1.40057},

>    {x -> 1.8183}}

In[3]:= Length[%]

Out[3]= 12
```

We should add a doctest for this, to check the bug does not come back once it gets fixed.

Dave

### comment:10 Changed 2 years ago by jdemeyer

• Milestone changed from sage-4.7 to sage-4.7.1

### comment:11 Changed 20 months ago by jdemeyer

• Priority changed from blocker to critical
• Milestone changed from sage-4.7.2 to sage-4.7.3

### comment:12 follow-up: ↓ 13 Changed 20 months ago by zimmerma

I'm concerned by the fact that this ticket was delayed twice (first from 4.7 to 4.7.1, now from 4.7.2 to 4.7.3, btw why is there no trace of the change from 4.7.1 to 4.7.2) and moreover the priority is now changed from blocker to critical, without any justification.

Is Carl Witty (who wrote the real_roots code as far as I know) still out there?

Paul

### comment:13 in reply to: ↑ 12 Changed 20 months ago by jdemeyer

I'm concerned by the fact that this ticket was delayed twice (first from 4.7 to 4.7.1, now from 4.7.2 to 4.7.3, btw why is there no trace of the change from 4.7.1 to 4.7.2) and moreover the priority is now changed from blocker to critical, without any justification.

Well, it seems nobody really cares much to solve this problem. I have glanced at the code and it's also hard to understand what is going on.

Personally I don't think this should a "blocker". Simple mathematical bugs are almost never blockers. This real_roots() code is also not the default real-root computing code. If you simply do pol.roots(ring=RR), this will use PARI and not this buggy code (there could still be bugs in PARI though).

Is Carl Witty (who wrote the real_roots code as far as I know) still out there?

AFAIK, he is not really into Sage development anymore.

### comment:14 Changed 20 months ago by zimmerma

Well, it seems nobody really cares much to solve this problem.

unless someone looks at it first, I will work on this ticket during the SageFlint? days in December.

Paul

### comment:15 Changed 19 months ago by jdemeyer

• Milestone sage-4.7.3 deleted

Milestone sage-4.7.3 deleted

### comment:16 Changed 12 months ago by dsm

• Milestone set to sage-5.1

Is it likely that anyone's going to look at this? If not, we should either put in a stopgap (now that we have them) or take a more drastic measure.

### comment:17 Changed 12 months ago by zimmerma

with the example in the description, over all seeds in [0,999], the number of roots is incorrect only for 1 and 304. Unfortunately 1 is the default seed!

Paul

### comment:18 Changed 12 months ago by zimmerma

• Keywords roots, sd40.5 added; roots removed
• Status changed from new to needs_review
• Authors set to Paul Zimmermann

I found out that reducing the Bernstein coefficients to integers fixes the two instances of this bug (the one in the description and the one from comment 8). Maybe it does not solve all instances, but it is better than nothing, and it should not hurt to work with integers instead of rationals.

I suggest reviewers try to find other instances, for example with:

```from sage.rings.polynomial.real_roots import real_roots
R.<x> = ZZ[]
B=10^11
for i in range(10^6):
f = R.random_element((i%13)+1,-B,B)
print f
sys.stdout.flush()
n = len(real_roots(f,seed=-1))
for s in range(1000):
if len(real_roots(f,seed=s))<>n:
print f, n, s
sys.stdout.flush()
raise ValueError
```

Paul

### comment:19 Changed 12 months ago by dsm

Given how sensitive the problem is I think we're going to need to know why whatever we find fixes the problem before we're confident that we've fixed it. Here I'm not sure that it does. Digging up some of my old notes, even with the patch, we have problems we shouldn't:

```sage: from sage.rings.polynomial.real_roots import real_roots
sage: R.<x> = ZZ[]
sage: f = 535944*x^12 - 3312356*x^11 + 8042902*x^10 - 9521317*x^9 + 5204589*x^8 - 461471*x^7 - 682572*x^6 + 174775*x^5 + 25797*x^4 - 5793*x^3 - 473*x^2 + 17*x + 1
sage: len(f.roots(RR))
12
sage: len(real_roots(f, seed=1))
12
sage: # but..
sage: f.roots(RR)
[(-0.279345282327294, 1), (-0.169895157384578, 1), (-0.0836442589128367, 1), (-0.0391184461806829, 1), (0.0525277889691767, 1), (0.217266273158562, 1), (0.608742206578326, 1), (0.743349076348040, 1), (0.955058967243006, 1), (0.956604779629507, 1), (1.40057197815310, 1), (1.81829644637645, 1)]
sage: [len(real_roots(f, seed=1, bounds=(-1, 7+1/2**k))) for k in [0..200]]
[12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 10, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12, 12, 12, 10, 12]
```

And we see a clearly periodic pattern. For me real_roots.pyx takes a long time to compile so I got frustrated chasing this clue, but changing the bounds like this -- outside the range of all of the roots! -- should have no effect. If we could figure out why this happens I think we'd have it.

### comment:20 Changed 12 months ago by zimmerma

I've worked a little more on this. It seems the _count_variations function is called with a wrong usign (sign in x=1 if I understand correctly) in the following case:

```enter _count_variations, self= <IBP: ((67412088, 41929949, 16504578, -2435470, -4640) + [0 .. 18683)) * 2^-10 over [3921/16384 .. 8017/32768]; usign 1; level 3; slope_err [-41131.637071203295 .. 41131.637071203295]>
```

and then returns 2 changes of sign whereas there is evidently only one in the coefficients.

Paul

### comment:21 Changed 12 months ago by zimmerma

I believe I've found the reason for the bug: the function down_degree decreases the degree of the (interval) polynomial, but keeps the value self.usign. With the example in the description and seed=1, down_degree is called with:

```enter down_degree, self= <IBP: ((34519150034, 25826190301, 17130808432, 9120071648, 2724478006, -830629355, 3694263) + [0 .. 21299)) * 2^-19 over [3921/16384 .. 8017/32768]; level 2; slope_err [-136.91992834614365 .. 136.91992834614365]> usign= 1
```

and the polynomial is reduced to:

```coeffs= (67412088, 41929949, 16504578, -2435470, -4640) n= 5
```

where it is clear that usign=1 is wrong, since the sign of -4640 is -1.

Paul

### comment:22 Changed 12 months ago by zimmerma

the new patch fixes the bug. Ready for review. Thanks Douglas for pushing me to really isolate the reason of the bug.

Paul

### comment:23 Changed 12 months ago by zimmerma

I just checked: the example from comment 19 passes too with the new patch.

Paul

### comment:24 Changed 12 months ago by dsm

Looks like nice work! I think we're going to have to add way more testing, though (I can do some of it as a reviewer patch this weekend.) I have a collection of failing cases and I want to make sure that more of them work; we should specifically doctest your diagnosis, not just the outcome; and I want to test that the bound stuff now works.

Now that we know the problem it should be much easier to construct more test cases, too.

### comment:25 Changed 8 months ago by jdemeyer

• Status changed from needs_review to positive_review
• Reviewers set to Jeroen Demeyer
• Milestone changed from sage-5.4 to sage-5.5

Since the reviewer patch didn't happen, let's put positive_review here. You can still add more examples on a new ticket.

### comment:26 Changed 7 months ago by jdemeyer

• Status changed from positive_review to closed
• Resolution set to fixed
• Merged in set to sage-5.5.beta0

### comment:27 Changed 7 months ago by jdemeyer

• Description modified (diff)
• Merged in changed from sage-5.5.beta0 to sage-5.4.1.rc0

### comment:28 Changed 7 months ago by jdemeyer

• Milestone changed from sage-5.5 to sage-5.4.1
Note: See TracTickets for help on using tickets.