Opened 11 years ago

Closed 10 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 sage-5.4.1.rc0 Paul Zimmermann Jeroen Demeyer N/A

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

### comment:1 follow-up: ↓ 2 Changed 11 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 11 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 11 years ago by jdemeyer

• Description modified (diff)

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

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

### comment:5 Changed 11 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 11 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 11 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 11 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 11 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 11 years ago by jdemeyer

• Milestone changed from sage-4.7 to sage-4.7.1

### comment:11 Changed 11 years ago by jdemeyer

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

### comment:12 follow-up: ↓ 13 Changed 11 years 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 11 years 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 11 years 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 11 years ago by jdemeyer

• Milestone sage-4.7.3 deleted

Milestone sage-4.7.3 deleted

### comment:16 Changed 10 years 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 10 years 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 10 years ago by zimmerma

• Authors set to Paul Zimmermann
• Status changed from new to needs_review

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 10 years 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 10 years 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 10 years 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 10 years 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 10 years ago by zimmerma

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

Paul

### comment:24 Changed 10 years 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 10 years ago by jdemeyer

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

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 10 years ago by jdemeyer

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

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

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

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

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