• March 01, 2021, 10:52:40 AM

Login with username, password and session length

Author Topic:  Round-off errors in polynomial evaluation  (Read 799 times)

0 Members and 1 Guest are viewing this topic.

Offline marcm200

  • 3c
  • ***
  • Posts: 925
Round-off errors in polynomial evaluation
« on: May 03, 2020, 09:21:52 AM »
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.

Code: [Select]
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

Code: [Select]
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
Code: [Select]
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:
Code: [Select]
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:

Code: [Select]
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/
« Last Edit: May 03, 2020, 10:50:09 AM by marcm200, Reason: clarification »


Offline marcm200

  • 3c
  • ***
  • Posts: 925
Re: Round-off errors in polynomial evaluation
« Reply #2 on: May 04, 2020, 01:19:24 PM »
@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).

Offline Adam Majewski

  • Fractal Frogurt
  • ******
  • Posts: 452
Re: Round-off errors in polynomial evaluation
« Reply #3 on: May 04, 2020, 05:14:25 PM »
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?

maybe slowly escaping points:
Checking one point if escapes in parabolic case :

for n = 34 take about one hour ( 5 595 seconds )
for n = 40 take about one day
for n = 45 take about one month
for n = 50 take about one year

https://en.wikibooks.org/wiki/Fractals/Mathematics/Numerical#Parabolic_case

Offline Adam Majewski

  • Fractal Frogurt
  • ******
  • Posts: 452
Re: Round-off errors in polynomial evaluation
« Reply #4 on: May 04, 2020, 05:21:16 PM »
It is probably numerical error. If yes , what type ?
Answer : Round-off ( as in the title ) (:-))


What about http://arblib.org/ ?
« Last Edit: May 04, 2020, 09:48:13 PM by Adam Majewski, Reason: answer »

Offline marcm200

  • 3c
  • ***
  • Posts: 925
Re: Round-off errors in polynomial evaluation
« Reply #5 on: May 05, 2020, 04:01:30 PM »
What about http://arblib.org/ ?

Yes, that would work. What made me stop for a bit when finding that example here, was, that usually when calculating with computer number types, using a higher precision result and comparing it to a lower precision one is a quick method to assess whether the result is (somewhat) correct. But seeing this now puts that method way down my list - and it makes me wonder, how many errors of that type are in point-sampled images or other results like the area of the Mandelbrot set.

Fortunately for me since I stumbled upon the Figueiredo algorithm last year (or more precisely, claude pointed me towards it), I'm mostly computing with interval arithmetics and dyadic fractions and use point-sampling only for artistic purposes.

I wonder how accurate Mandelbrot's first image of his set was. Does anyone know, which type of computer he used and what number type?


Offline claude

  • 3f
  • ******
  • Posts: 1784
    • mathr.co.uk

Offline marcm200

  • 3c
  • ***
  • Posts: 925
Re: Round-off errors in polynomial evaluation
« Reply #7 on: May 06, 2020, 05:43:30 PM »
EDIT 2020 05 07: The 5-bit part of the image below is not correct (please see post 10):

What about http://arblib.org/ ?
Do arbitrary precision libraries deal with periodic numbers like 1/3 in a specific way - or are they assigned as the binary equivalent to 0.3333... with the maximum number of digits?

@claude: Interesting link. It reads almost like a detective story, trying to catch the cukprit (the grid used for the complex plane tiling).

I tried here the opposite to arbitrary precision: lower precision and simulated binary numbers with fewer than 52 bits fractional. Visually the Mset starts to look different at around 10 bit precision, but file-wise, 50 bits precision has pixels differently colored than at 52 bits. So on one side it's human-visually stable, on the other side it's unstable just using a couple of fewer bits.

Technical details:
  • depicted is the 2-square, image computed by point sampling, 100 maxit
  • colored by log(escape iteration)/log(1000) as index to a heat-map
  • lower precision of a double-precision value B was simulated as: lowprec(B) := truncate(B*2^bits) / 2^bits
  • evaluation works in a tree (a,b,c low prec numbers): (a op1 b) op2 c := lowprec(a op1 lowprec(b op2 c))

« Last Edit: May 07, 2020, 09:53:29 AM by marcm200 »


Offline hobold

  • Fractal Fluff
  • *****
  • Posts: 397
Re: Round-off errors in polynomial evaluation
« Reply #9 on: May 06, 2020, 09:20:32 PM »
Do arbitrary precision libraries deal with periodic numbers like 1/3 in a specific way - or are they assigned as the binary equivalent to 0.3333... with the maximum number of digits?
Most libraries use a (much) larger number of bits and do the usual arithmetic on them, including round-off errors. But there are a few specialized libraries that work with ratios of (arbitrarily large) integers. I haven't seen the latter used for numerical bulk computation, though; I encountered them in whitepapers on computation with mathematical guarantees.

Offline marcm200

  • 3c
  • ***
  • Posts: 925
Re: Round-off errors in polynomial evaluation
« Reply #10 on: May 07, 2020, 09:51:16 AM »
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?


Offline Adam Majewski

  • Fractal Frogurt
  • ******
  • Posts: 452

Offline gerrit

  • 3f
  • ******
  • Posts: 2333
Re: Round-off errors in polynomial evaluation
« Reply #12 on: May 07, 2020, 04:51:00 PM »
@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.

Offline claude

  • 3f
  • ******
  • Posts: 1784
    • mathr.co.uk
Re: Round-off errors in polynomial evaluation
« Reply #13 on: May 07, 2020, 05:49:47 PM »
But that's not the same as saying you can't make a reliable picture with pixels.
https://doi.org/10.1002/malq.200310124 says it depends on the "famous hyperbolicity conjecture"

Offline gerrit

  • 3f
  • ******
  • Posts: 2333
Re: Round-off errors in polynomial evaluation
« Reply #14 on: May 07, 2020, 08:59:02 PM »
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.


clip
A Round Thing

Started by Bill Snowzell on Fractal Image Gallery

2 Replies
217 Views
Last post March 09, 2018, 09:55:20 AM
by Bill Snowzell
xx
CPU: order of evaluation of an arithmetic expression

Started by marcm200 on Programming

5 Replies
299 Views
Last post February 12, 2020, 09:35:38 AM
by hobold
xx
Several Problems / Errors

Started by PMH on Kalles Fraktaler

5 Replies
470 Views
Last post December 07, 2017, 04:43:37 PM
by PMH
xx
Errors when trying to upload video

Started by LionHeart on Forum Help And Support

6 Replies
240 Views
Last post April 13, 2020, 02:36:08 AM
by LionHeart
xx
compile errors with release 2-2.13.

Started by piotrv on Mandelbulber

3 Replies
359 Views
Last post April 01, 2018, 07:42:03 AM
by buddhi