## Bug #1171

### RealRoots: first point is sometimes wrong?

**Description**

It seems that sometimes (very rarely) the first root given by `RealRoots`

is wrong

/**/ use QQ[x]; /**/ f := x^25 +1500*x^24 -3472*x^23 +75*x^22 -2888*x^21 +12984*x^20 +(-22927/4)*x^19 +11578*x^18 -22029*x^17 +(18981/2)*x^16 -21819*x^15 +29021*x^14 -12855*x^13 +25232*x^12 -25419*x^11 +4494*x^10 -20702*x^9 +14894*x^8 -1920*x^7 +12502*x^6 -92*x^5 +1700*x^4 -6628*x^3 -616*x^2 -520*x +920; /**/ RR := RealRoots(f, 1/1000); /**/ DecimalStr(subst(f, x, RR[1].inf)); -6710019372583766710138791136389510647991731589224733661634966329959500813.670 /**/ DecimalStr(subst(f, x, RR[2].inf)); 0.034 /**/ DecimalStr(subst(f, x, RR[3].inf)); -0.418

**Related issues**

### History

#### #1 Updated by Anna Maria Bigatti 4 months ago

I found this bug while debugging `ApproxSolve`

.

I had noticed that, when the answer was wrong, it was just on the first point.

#### #2 Updated by Anna Maria Bigatti 4 months ago

It seems that actually the problem is in the instability of the roots: the approximate root is probably close to the real solution, but the evaluation is not close to 0.

/**/ f := x^25 +1500*x^24 -3472*x^23 +75*x^22 -2888*x^21 +12984*x^20 +(-22927/4)*x^19 +11578*x^18 -22029*x^17 +(18981/2)*x^16 -21819*x^15 +29021*x^14 -12855*x^13 +25232*x^12 -25419*x^11 +4494*x^10 -20702*x^9 +14894*x^8 -1920*x^7 +12502*x^6 -92*x^5 +1700*x^4 -6628*x^3 -616*x^2 -520*x +920; /**/ [DecimalStr(subst(f, x, rr.inf)) | rr in RealRoots(f, 10^(-50))]; ["-1.095", "0.000", "-0.000", "0.000", "-0.000"] /**/ [DecimalStr(subst(f, x, rr.inf)) | rr in RealRoots(f, 10^(-100))]; ["-0.000", "0.000", "-0.000", "0.000", "-0.000"]

Another example on this pathology.

/**/ f := x^25 -3472*x^23 +11578*x^18 +x^7 +92*x^5 +x^4 -x^3 -x^2 -x +920; /**/ [DecimalStr(subst(f, x, rr.inf)) | rr in RealRoots(f, 10^(-100))]; --> 0 0 0 /**/ [DecimalStr(subst(f, x, rr.inf)) | rr in RealRoots(f, 10^(-10))]; --> ["-1749714807099581247924392.282", "0.000", "-1882915635399905369900820.013"]

#### #3 Updated by John Abbott 4 months ago

**Status**changed from*New*to*In Progress***% Done**changed from*0*to*50*

One can estimate the "stability" of the root by computing the value of the derivative at that root. By Taylor expansion we have f(alpha+eps) = f(alpha) + eps*deriv(f)(alpha) + other terms.

Probably we would want an upper bound for abs(Df(beta)) where beta is a "good approximation" to alpha.

#### #4 Updated by John Abbott 4 months ago

**Related to***Feature #1173: Upper bound for value of poly in an interval*added

#### #5 Updated by John Abbott 4 months ago

It is not as simple as I first thought. While the first order approximation I mentioned in comment 3 is correct, it is only a first order approximation.

I made an implementation, and I tried the example in comment 2. That worked OK.

But with the polynomial **f^2+f** problems arose, presumably because the first order approximation was not good.

It would be nice to have a fn which returns intervals with the guarantee that **max(abs(f(x)) | x in [a,b]) < epsilon** where **epsilon** is the bound specified by a user.

#### #6 Updated by John Abbott 4 months ago

- call
**RealRootsSmallValue(f, eps)** - result is a list of closed intervals whose interiors are disjoint such that

- width of interval is less than**eps**(less than or equal?)

- value of**f**over each interval is less than**eps**(less than or equal?)

It is permitted to return intervals which are "too small".

What do you think? Suggestions for the fn name?

#### #7 Updated by John Abbott 2 months ago

I am not so happy with the suggested interface in comment 6: the parameter `eps`

is doing double duty (a bound for the accuracy of the root, and a bound for the value of the function). Having two separate parameters for these bounds would be semantically neater (but possibly more awkward for the caller?)

Comments? Suggestions?