Hacker News new | past | comments | ask | show | jobs | submit login
How Mandelbrot set images are affected by floating point precision (github.com/profjski)
271 points by todsacerdoti 9 months ago | hide | past | favorite | 95 comments



From a cursory look at the code, I don't see any use of fused-multiply-add (FMA) which would likely help with precision issues in a number of places with the float version. The "problem in more detail"[1] readme specifically calls out computation of `x^2 - y^2` as a source for error, and that has methods that dramatically reduce error with FMAs[2].

[1] https://github.com/ProfJski/FloatCompMandelbrot/blob/master/... [2] https://pharr.org/matt/blog/2019/11/03/difference-of-floats


The author needs to learn about techniques for rendering deep zooms into the Mandelbrot set. Since 2013 it has been possible to render images that are 2^-hundreds across using mostly double precision arithmetic, apart from a few anchor points calculated with multiprecision (hundreds of bits) arithmetic.

The deep zoom mathematics includes techniques for introspecting the iterations to detect glitches, which need extra multiprecision anchor points to be calculated.

https://mathr.co.uk/blog/2021-05-14_deep_zoom_theory_and_pra...

https://dirkwhoffmann.github.io/DeepDrill/docs/Theory/Mandel...

In my experience of non-deep-zoom rendering, and contrary to the authors arguments, period detection works well for speeding up renders. It appeared to be fairly safe from false positives. https://dotat.at/@/2010-11-16-interior-iteration-with-less-p... https://dotat.at/prog/mandelbrot/cyclic.png


I've always wondered how those "Ten minute Mandelbrot zoom" videos worked, because there's no way double would last that long at a tolerable zoom rate.

The perturbation technique is interesting. Calculating just a few points with super high precision and then filling the pixels in between by adding an offset and continuing with lower precision halfway though the calculation seems plausible at a glance, but I'll have to read that more carefully later.

Thanks for sharing!


I've had some great success using posits for fractals. Unfortunately, soft(ware) posits so rather slow. But mathematically the results were great.


super interesting thanks, really.


But the goal of the study is expressly to highlight precision issues. Using such techniques wouldn't help to make them surface.


I wish someone made up a better explicit FMA syntax than std::fma(-c, d, cd) ... maybe something along the lines of ((-c * d + cd)) with special brackets or (-c) ⊠ d + cd with a special multiplication symbol.

And if only we could effectively use FMAs from javascript...


Mandelbrot set calculation is a curious intersection of math and computer science. I went down that rabbit hole once :)

Naive implementation is easy, but there are many ways to improve performance.


Strange thing to say since CS is a sub-field of math.


I wouldn't really say that's a universally held notion, especially in the modern world. Historically, yes, but these days it stands more on its own. In similar vein, you could say biology is a subfield of physics, but most people don't think of them that way.


To a librarian, Computer Science and Library Science are just different aspects of information systems, which are (a librarian need hardly tell HN) humanity's highest calling, thus justifying their shared Dewey classification of 0.


Interesting! And "Science" is all the way down in Class 500.

https://en.wikipedia.org/wiki/List_of_Dewey_Decimal_classes


That's so nice, it brings a tear to the eye. I hope we live up to librarians' expectations.


It's not really philosophically, though. As a pure mathematician, I would not consider CS a subfield of math. The spirit of math is really studying abstract axiomatic systems, and while CS in theory is sort of like that, in actual practice, CS is too tied and focused on the actual instantiation of computers to really have the spirit of math.


Many would disagree with your characterisation of pure mathematics, computer science, or both. And I am certain that computer science yields nothing to pure mathematics in her pursuit of abstraction.


All I know is that when I read computer science, it doesn't feel like math at all. Other people can characterize it however they choose.


Most CS in practice is really software engineering


I think that statement held true before widespread availability of computing; when most computer science was really theoretical. I once skimmed a few chapters of a graduate-level textbook on Category Theory and realized that it was the foundation of object-oriented programming.

The biggest issue is that a lot of "Computer Science" is really applied software engineering; much like confusing physics and mechanical engineering.

Or, a different way to say it: Most students studying "Computer Science" really should be studying "Software Engineering."

More practically, I have a degree in Computer Science. When I was in school, it was clear that most of my professors couldn't program their way out of a paper bag, nor did they understand how to structure a large software program. It was the blind leading the blind.


Meanwhile at my ABET acredited, 2nd rate state school, my Computer Science college professors were talented people who had been doing this shit for decades and clearly understood what they were talking about.

I had ONE class that was about Software engineering.

Meanwhile I had an entire years worth of curriculum that was just "Go take various unrelated science and math classes so you have a strong understanding of the fundamentals in both science and math"

People so often generalize their very specific college experience to the entire world. Meanwhile you'd be lucky to find a consistent college experience just from crossing state lines.


All fields are subfields of polymathy.


I thought you're joking, but nope.

https://en.wikipedia.org/wiki/Polymath


Of course these terms aren’t well-defined, but some (me) would say both math and CS are subfields of “formal systems” (but “formal systems” could be named CS, which would put it at the top of the hierarchy)


The Joy computer language and Scheme with SICP are prime examples.

Reading SICP it's the 'easy way' and the Joy docs plus trying to code even simple as a quadratic eqn solver it's a deep, hard task which requires lots of knowledge on cathegory theory.


Software engineering has as much to do with sociology as it does with math.

Systems software has to deal with a lot of physics and engineering stuff (speed of light, power, heat, mean time to component failure, etc, etc.)


Software Engineering and Computer Science don't necessarily mean the same thing (even though at this point, most schools with CS programs just teach software engineering anyway).


…and I thought math branched off from computing science sometime in the 19th century. Before that, what was called „math“ was mostly algorithms to compute stuff.


This post prompted me to look into the theoretical question of the computability of the Mandelbrot set. And I found this excellent answer: https://math.stackexchange.com/q/4611850


Another way to explore floating point representation can also be to compare rounding modes, as this can impact a lot of what is generated for floating point operations. Most systems are round to nearest even, but round to zero, round towards positive and round towards negative also exist. Bonus points if you also check results for x87 floats.


For Mandelbrot, presumably you can do exact rational arithmetic (all it needs is multiplications and add, and you presumably start with a floating-point coordinate, which is eminently rational). It will be very slow with many iterations, of course, since the fractions will rapidly spin out of control.

Edit: I see the supplementary material addresses this: “One could attempt to avoid the problem by using arbitrary-precision floating point representation, or rational representation (using two integers, for numerator and denominator, i.e., keep everything as fractions), but that quickly leads to so many significant digits to multiply in every iteration that calculations become enormously slow. For many points in a Mandelbrot Set Image, the number of significant digits in a trajectory tends to double with each iteration, leading to O(n^2) time complexity for evaluating one point with n iterations. Iterations must then be kept very low (e.g., 20-30 on a modern PC!), which means that many higher-iteration points cannot be evaluated. So accumulated round-off error is hard to avoid.” 20-30 sounds crazy low even for O(n²), though.

Edit 2: This part seems wrong: “Cycle detection would work if floating point quantities could be compared exactly, but they cannot.” They certainly can. If you detect a cycle, you will never escape _given your working precision_. Floats can cycle just as well as fixed-point numbers can. But I suppose true cycles (those that would be cycles even in rationals) are extremely unlikely.


“One could” indeed:

https://www.fractint.org/

From wikipedia:

> The name is a portmanteau of fractal and integer, since the first versions of Fractint used only integer arithmetic (also known as fixed-point arithmetic), for faster rendering on computers without math coprocessors. Since then, floating-point arithmetic and arbitrary-precision arithmetic modes have been added.


I spent a lot of time with fractint in the late 1990s. Great memories!


Fractint is fixed-point, not exact rationals.


Hm, what about intervals of rationals, and then approximating the intervals with a slightly wider interval in order to make the denominators smaller? I guess the intervals for the real and imaginary components might grow too quickly?

if 0 < a_L < a < a_R , 0 < b_L < b < b_R , then (a + b i)^2 = (a^2 - b^2) + 2 a b i , and a_L^2 - b_R^2 < a^2 - b^2 < a_R^2 - b_L^2 and 2 a_L b_L < 2 a b < 2 a_R b_R ,

uh,

(a_R^2 - b_L^2) - (a_L^2 - b_R^2) = (a_R^2 - a_L^2) + (b_R^2 - b_L^2) = 2 (a_R - a_L) ((a_L + a_R)/2) + 2 (b_R - b_L) ((b_L + b_R)/2)

And, a_R b_R - a_L b_L = a_R b_R - a_R b_L + a_R b_L - a_L b_L = a_R (b_R - b_L) + (a_R - a_L) b_L

So, looks like the size of the intervals are like, (size of the actual coordinate) * (size of interval in previous step), times like, 2?


"20-30 sounds crazy low even for O(n²), though."

I mean, 30 means 2^30 sig digits. That means any operation on those are a bit expensive in pure memory throughput. Squaring that, let's use an FFT then it's just O(n log n) for the next step, with n being number of digits here. Or 30 * 2^30 ops, let's call it 2^35 ops. We're at 32 MFlops for a single dot. Let's only do a 1280*768 image, that's 32 GFlops to compute that image.

That's a perfectly nice rate, let's pretend it's half a second because all memory access was just perfectly sequential. (It's not, it's more).

Each additional step more than doubles the time spent. Which means even at 6 more iterations, we're staring at half a minute compute time. Nobody waits half a minute for anything. (Unless it's a web page downloading JS)

Now add the fact that I lied, memory access can't be fully sequential here. But even if it was, we're at least spilling into L3 cache, and that means some 40 cycles penalty every time we do. We'll do that a lot. So multiply time with 10 or so (I'm still hopelessly optimistic we have some efficiencies.), and suddenly even 30 iterations cost you 5 seconds.

And so, yes, it's really expensive even for modern PCs.

If you have a GPU, you can shuffle it off to the GPU, they're pretty good at FFTs - but even that's going to buy you, if we're really really optimistic, another 20 or so iterations before it bogs down. Doubling digits every cycle will make any hardware melt quickly.


Ah, of course, not only does it require O(n²) work to multiply, but n grows to n² for each iteration (n bits -> 2n bits). So that's why you get O(2^n) _total_ time.


afaik all ieee floating point systems allow configuring rounding mode, for example with fesetround or equivalent https://pubs.opengroup.org/onlinepubs/9699919799/basedefs/fe...


I had never realized the necessity of nearest even.

I thought that 0.5 ought to be rounded to 1, mirroring 0.0 rounded to 0 !

But I can now see how this only works for reals, while anything non-continuous (even rationals ?) could end up with compounding errors...


I've heard it described as bank-rounding


In Hamiltonian Chaos you draw Poincaré sections of systems like this one

https://mathworld.wolfram.com/Henon-HeilesEquation.html

and see these chaos zones that look like a field of random dots that probably miss a lot of the structure. The reason why we know those areas are dense with unstable periodic orbits is because they are also dense with periodic orbits. It doesn’t seem feasible to accurately track a system for long enough to believe those.


The second[1] C program I read in my life was Fractint[2]. It is called that way because it calculated fractals with integer math[3].

At that time integer math was not only faster than floating point, but the only thing available to most of us. Processors like the 80386 did not come with any floating point facilities whatsoever. The FPU was a separate unit (aptly called a coprocessor at the time) which was expensive and hard to obtain where I lived.

[1] The first one was DKBTrace, which lives on as POVRay. Thank you David K. Buck, where ever you are now. I learned a lot from you.

[2] https://fractint.org/

[3] Technically we should call it fixed point I guess, but in the end it was all done with the integer instructions of the CPU. So integer math is not completely wrong.


I still have the entire manual that I printed and bound circa 1996. I saw a program on PBS about fractals and we had just got dial up internet. I found fractint and cut my programming teeth with that and a copy of Borland that I convinced my parents to buy for me. Running on a 75MHz Compaq. Being a teenager at that time was just magical.

Edit: Just noticed you mentioned Povray too. Are you me? Wow, lol.


Haha, I also still have the manual which I printed on my Epson LQ-400 24 needle printer when no one was home because it was so loud and took forever.


I think he's in Ontario, but you can ask him yourself: david at simberon dot com ;)


"Processors like the 80386 did not come with any floating point facilities whatsoever."

The 80286 needed a 80287 - I bought one to run AutoCAD (1MB RAM!)

I can't quite remember the exact order of events but either or both 80386 and 80486 had a co-pro initially. Then Intel re-invented the IBM magic screwdriver. They dropped the 80 prefix and added SX or DX and an i prefix for a laugh. SX had no co-pro or it was disabled and DX was the full fat model.

They really started to screw the punter with the 586 ... Pentium. The trend has continued - Core iBollocks.


They offered the 80387 co-processor a few years after the 80386 came out

ref https://en.wikipedia.org/wiki/X87#80387


I had an 80386DX (or maybe i386DX, but DX for sure) and it definitely had no FPU. Always wanted one so badly.


I was going to correct you, and say the DX had a coprocessor, but the SX did not.

Apparently, that’s incorrect, and the difference is bus width.

However, at the end of its life, my 386 motherboard was running one of these at 100MHz, with a coprocessor on the daughterboard in the main CPU socket:

https://en.m.wikipedia.org/wiki/Evergreen_Technologies

(I think I had an IBM 486, not a Cyrix, fwiw.)

Sadly, I’d moved to Linux at that point, and the processor needed a special MSDOS binary to flip into 100MHz mode after boot.

Even with the handicap, it booted into X11 + enlightenment at 25MHz (or maybe 33MHz) faster than it could boot into Windows 3.11 at 100MHz.


PLCC i80387 or compatible can be had from eBay for less than $20, and PGA generally for less than $30 USD.


For that price you can just get a 486 DX.


...which won't fit in an 80387DX socket...


If I remember correctly the 486 was just basically, at first, the 386 with the FPU built in. There may have been a couple of new instructions? On that last part I am a bit fuzzy.


The 8086/8088/80186/80188, 80286 & 80386 all required an external chip for floating point.

The 80386DX was the full 32-bit data bus version. The 80386SX was a cost reduced version with a 16-bit data bus (akin to the difference between the 8086 & 8088).

The 80486DX had an FPU; the 80486SX had a disabled or missing FPU (depending on vintage). You could pair the 80486SX with an 80487 FPU, which was just a full 80486DX in a different pinout that took over from the crippled part.

While some Intel doco did start calling them Intel386 & Intel486, they didn't drop the 80 until Pentium. Intel lost a court case trying to prevent clone makers from calling their chips '80486' or whatever because the judge ruled you can't keep your competitor from using a number.

Weitek made floating point units that paired to the 386 and 486. Some motherboards supported this, and a small amount of software did (early Autocad).

Correcting one down thread, the 486 is more than just a 386+387. There a cache, it's pipelined, and there are some changes to the ISA. 4x the transistors.


Yeah, Fractint was great. The VGA color palette cycling was trippy.

In the mid 1980's I typed in a Mandlebrot set program from a magazine. It was for the Apple 2, written in BASIC. It took around 12 hours to generate a 280x192 monochrome Mandlebrot.


Wow! I used to use fractint ages ago, and I thought "tint" was for colors...


If you need to use 128-bit floats for something, note that your processor probably supports them natively and that there are standard C/C++ types for them.

(The boost thing the article links says their quad precision type is “functionally equivalent” to standard types, but uses a different bit layout.)


Your processor probably doesn’t support quad precision floats in hardware unless it’s an IBM mainframe. https://en.m.wikipedia.org/wiki/Quadruple-precision_floating...


Power also supports it, and https://www.raptorcs.com/ aren't mainframes.


Which prompts the question, what is quad precision floating point arithmetic needed for that IBM considered it prudent to implement it in hardware of a mainframe? Who are the customers paying for such? Naively perhaps, I would expect that financial calculations (ought to) rely on fixed-point arithmetic and afaiu simulations of physical systems get away with double precision, if one is careful.


>afaiu simulations of physical systems get away with double precision, if one is careful.

Doesn't this depend very much on the physical system in question. I could imagine that simulating chaotic systems needs a lot of precision.

Simulation of n-body problems and others are mentioned [1].

[1] https://www-zeuthen.desy.de/acat05/talks/de_Dinechin.Florent...


I've read this a few places recently but I can't figure out how? https://godbolt.org/z/Mbjs73qar -- using doubles I get AVX instructions but using 128-bit floats calls some very hefty runtime support functions.


There is no true quad support on x86. GCC implements double-double quads[1] via software emulation on many targets.

[1] https://en.wikipedia.org/wiki/Quadruple-precision_floating-p...


Standard C/C++ types may still use a software implementation behind the scene, and this is almost always the case for std::float128_t.


A related question I've wondered is, can you rearrange the formulae to reduce precision loss, without having to actually use more bits of precision in the data type?


Long time since I looked at this but you might be interested in:

https://mathr.co.uk/blog/2021-05-14_deep_zoom_theory_and_pra...

IIRC it is calculating a few points at high precision and then estimating the rest.


Mandelbrot specifically, or computations in general? For the latter, absolutely yes. See Herbie for automated tool for doing it https://herbie.uwplse.org/


The former


Did anyone else think it was interesting that the difference between the two floating point images looked like a Julia set? It’s probably trivial, but I liked that.


It's essentially just the boundary of the Mandelbrot set which makes sense as that's where the calculation is the most unstable.


Mandelbrot and Julia sets are related:

https://www.karlsims.com/julia.html


The popular Mandelbrot zoom-in videos all seem to develop rotational symmetry before long, and the order of that symmetry only increases over time. I've assumed that this is artifact of precision errors compounding.


It’s because they usually zoom in on a high-order hyperbolic component https://en.wikipedia.org/wiki/Mandelbrot_set#Hyperbolic_comp...


https://www.youtube.com/watch?v=aP0Y1uAA-2Y

A good lecture by John L. Gustafson on beyond floating points?


25 years ago I did Mandelbrot with fixed-point integer, and I remember zooming in enough and eventually hit the end of the fractal due to limited resolution.


I fondly remember Fractint (https://www.fractint.org/). And these days I'm amazed at what it could calculate with the very limited CPU cycles at that time. I'm not sure how to compare a 386 running at 25MHz with my AMD 5950x that I have now.


Fractint did a lot of clever tricks, especially for the Mandelbrot fractal. For example, I remember they exploited that the Mandelbrot set is connected. So that means that if you found that the border of any rectangle is either completely inside or completely outside the set, you don't have to check the inside.


Irrelevant nostalgic note: I implemented the Mandelbrot on the ZX81. Took just a quarter hour or so for the full glorious 64x48 pixels... :-)


Computing the color of a pixel is computing the (averaged) integral of some iteration counting function. It looks like there are areas where this function is pretty chaotic -- isn't the problem very ill-conditioned? If so adding more precision may not even help?


Seems like this algorithm does not do that ?

https://github.com/ProfJski/FloatCompMandelbrot/blob/master/...

> [...] Instead, we put pixels in a one-to-one relation to the points they represent. You can think of this as mapping some imaginary point on the pixel, such as its center or corner, onto the complex plane, and making the whole pixel represent that value.


I don't get the point of the project. So they show that rounding errors are amplified because a problem is ill-conditioned. Yeah, of course. The higher precision computation they do is probably equally wrong then.

Averaging / smoothing can sometimes help to make a problem better conditioned, so I expected they'd do something like that to get an accurate solution regardless.


It is equally instructive to compute the set on CPU and on GPU.


This reminds me of the issues with calculating the zeros of the Riemann zeta equation - your code design quickly becomes dominated by numerical implementation considerations


Some people, when they have a problem, think, “I know, I’ll use floating point.” Now they have 1.9999998 problems.


Funny story from today while briefly on hold for a support rep. Got the standard 'we'll call you back if you want blah blah' followed by - I kid not - "you're estimated wait time for the next representative is 0.58333333 minutes"

At the end of the call I suggested if she can log a problem for IT to include this and that they should either round up or use seconds (I suspect it was 35 seconds)


We had this issue with time estimates at work. JIRA field was minutes, which got printed as fractional hours elsewhere. Not good if you put 20 minutes. So suggestion was to use 15 minute intervals, but I pointed out that as long as you used multiples of 3 you'd get 2 digit results, so I'd round to multiples of 3 instead


Those are ridiculously fine grained tasks (or very precise estimations)?

Anything that takes less than an eg half an hour, I would probably just do, instead of making a ticket and an estimation. But I guess some companies are more bureaucratic than others?


Some numbers can be precisely represented using IEEE-754 floats. Two (2) is one of them.


Yes, but not every calculation that yields 2, also yields 2 when performed in IEEE-754 floats.


Yes, but any calculation that you'd use to count whole units things, like reasons, that would result in 2 would yield precisely 2.


Yes. The quip was only a silly joke.


I bet you have been waiting years to pull that one out of your pocket.

Well played sir! Nice shot man! :D


The document ("What Every Computer Scientist Should Know About Floating-Point Arithmetic") linked from the article was IMO even more interesting than the article: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.h...



Ironic that the last two letters in the tittle are truncated, "Mandelbr".

Must be a floating point rounding error.


And thats why I only use `fractint`! ;)


Code from "The Computational Beauty of Nature" managed to crash CWM under GNU/Linux+X. No issues with FVWM. WIth LibdieHard I managed to avoid the crash.

Maybe some ANSI C code from the 90's doesn't like current layouts...

On FP precision, sh has almost none, AWK often has up to 20 bit and 53 with gawk and the math library. Ditto with bc -l.

Calc has an arbitrary decimal precision and it's really fast. https://github.com/lcn2/calc




Guidelines | FAQ | Lists | API | Security | Legal | Apply to YC | Contact

Search: