I recently came across a polynomial that was new to me and that illustrates nicely what can go wrong with floating point arithmetics in computers.

Global optimization to prescribed accuracy
Ramon E Moore
Computers Math. Applic. Vol 21, 6/7 25-39, 1991
Chapter 3: Round-off errors
fulltext: https://www.sciencedirect.com/science/article/pii/089812219190158Z

This polynomial (EDIT: in variable a, if b is kept constant) evaluates very badly at some points - and all input numbers are accurately representable by double precision, so any error occurs due to round-off.

\[
333.75\cdot b^6+a^2\cdot (11a^2\cdot b^2-b^6-121b^4-2)+5.5b^8+\frac{a}{2b}
 \]

where a=77617, b=33096

checking double/float128 gives the same result

double:   f=1.17260394005317869492444060597
float128: f=1.17260394005317869492444060597

but that is not even correct in the sign. The correct answer is: -0.827396...

Interestingly, even WolframAlpha fails here with the standard settings
input  : 333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/(2b) where a=77617, b=33096
oputput: 1.18059×10^21

Only, when I force it to use rationals, it gives the correct output:
input : (333+3/4)*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+(11/2)*b^8+a/(2b) where a=77617, b=33096
output: -54767/66192 = -0.827396...

And even more interesting is, that if I change the b-value by 1 unit to 33095 - now all computations are agreeing (to some extent) and correct:

double f=-4.78339168666786449356823778558e+032
float128 f=-4.78339168666055425065308999647e+032
WA (fraction) -63322539148012414193286707611938758031/132380
f=-4.783391686660554025780836048643205773606284937301707 × 10^32

I hope I can find some points in polynomials for which I compute (non-IA) Julia sets that show the wrong escape behaviour.

Does anyone know of lower degree polynomials with such a behaviour?



Linkback: https://fractalforums.org/programming/11/round-off-errors-in-polynomial-evaluation/3460/

@Adam: Thanks for the links. It seems you have a large collection of references for a lot of topics. (I started one for myself a while ago for interesting things).

It is probably numerical error. If yes , what type ?
Answer : Round-off ( as in the title ) (:-))


What about http://arblib.org/ ?

https://math.stackexchange.com/questions/2358889/grid-spacing-iterations-used-in-the-1978-first-published-rendering-of-the-mande/3510304#3510304 is related, but does not have the answer to your specific question

Recall https://fractalforums.org/fractal-mathematics-and-new-theories/28/perturbation-theory/487/msg2358#msg2358

CORRECTION: The 5-bit-precision image from yesterday is not correct. I accidently commented out the lowprec-ing of the seed value and the image scaling factor. 5-bit would result in an empty image using a 2^10 pixel width, as scaling would be zero (4/2^10 = 2^-8 with 5 fractional bits) and the seed always -2-2i. Sorry.

@gerrit:

I read your post about the shadowing lemma, maybe you can help me with some questions:

Starting with a seed c and computing the orbit with round-off errors, the lemma states, there is a c2 close by that has this exact orbit (or an orbit arbitrarily close to it, that's not clear to me) using infinite precision. But one iteration later, it might be a different c2, right? And c and c2 do not necessarily have to have the same true fate, or do they?

So the lemma cannot help me about where my original c actually is going to, I could still have falsely classified it. Or am I missing something here?

@gerrit:

I read your post about the shadowing lemma, maybe you can help me with some questions:

Starting with a seed c and computing the orbit with round-off errors, the lemma states, there is a c2 close by that has this exact orbit (or an orbit arbitrarily close to it, that's not clear to me) using infinite precision. But one iteration later, it might be a different c2, right? And c and c2 do not necessarily have to have the same true fate, or do they?

Yes to all, specifically c and c2 are in the same pixel square.

So the lemma cannot help me about where my original c actually is going to, I could still have falsely classified it. Or am I missing something here?

There is no algorithm that can classify an arbitrary c (given with infinite precision). But that's not the same as saying you can't make a reliable picture with pixels.

Penrose has a lot to say about the issue of computability of the M-set in "The Emperor's New Mind: Concerning Computers, Minds, and The Laws of Physics (1989)". It seems the outside of the M-set is recursively enumerable but not recursive, which means you can verify algorithmically that a point that is outside is in fact outside, but you can't do that for any point inside. This uncomputability has not been proven however. And it's assumed that c itself is computable of course otherwise it's trivial.