5 ms·

Optimization Example: Mandelbrot Set (part 1)

74 points by RoyWard 3y ago

LorenDB 3y ago
Could you drop a link on your site to https://www.orange-kiwi.com/index.xml https://www.orange-kiwi.com/index.xml? Hugo generates that RSS feed by default, but you haven't linked it anywhere. Just a friendly heads up :)
- tsukurimashou 3y ago
  <link rel="alternate" type="application/rss+xml" title="Orange Kiwi" href="https://www.orange-kiwi.com/index.xml" /> yes just add this inside your <head>
deleted 3y ago
[deleted]
dist1ll 3y ago
Cool article. I think what the author calls "sheep race" is known as loop rotation. This allows you to mask off latency, at the cost of wasted computation at the tail of the loop.
- RoyWard 3y ago
  Do you have a reference for that? I might be looking in different places, but my reading suggests that loop rotation is where a conditional is moved to the end to avoid an unconditional branch, sometimes by the compiler. "Sheep race" with it's wasted computation is a little different - it's not a transformation it would be legal for a compiler to make, and possibly doesn't have enough use cases to really have its own well known name. Oddly enough, I've noticed that the compiler generated code will sometimes not do loop rotation, and and will have a conditional branch in the middle and unconditional at the end. I suspect this is because (a) unconditional branches are often effectively free on modern processors (if all the instructions are cached appropriately), and (b) there are all sorts of subtle instruction layout and ordering that a good compiler does to make things map efficiently to the underlying microarchitecture that means that non-intuitive orderings sometimes have a better cost function. I did find out that what I called "interleaving" is more properly called "software pipelining".
selimthegrim 3y ago
Are the streams problems good concurrency practice?
crazygringo 3y ago
Wow, this article opened my eyes to what is possible in optimizing code for a CPU. I was expecting this article to rehash the usual algorithmic tricks for speeding up Mandelbrot set calculation -- filling in boxes where the four corners already have the same color, checking for repeating iteration cycles to fill in black, etc. But I was very surprised to learn about how important CPU pipelines are, and therefore about "sheep race optimization" and interleaving, and then about the possibility of using vectorization on CPU. And in the end I'm utterly astonished that adding everything up results in the algorithm being calculated eight times faster. I would never have dreamed that such a level of speed improvement would be possible on such a simple, straightforward mathematical algorithm. Of course, in the real world everybody today is going to calculate the Mandelbrot set using the GPU (as the author notes -- this is just an exercise). But now it has me tremendously curious as to whether there are similar optimizations that can be applied on the GPU side? And especially whether those are accessible if done in WebGL/WebGPU? There are quite a number of Mandelbrot set generators available in the browser that use the GPU for calculations (e.g. [1]). [1] http://mandelset.ru/ http://mandelset.ru/ - https://github.com/iskandarov-egor/mandelset https://github.com/iskandarov-egor/mandelset
- QuadmasterXLII 3y ago
  "Are there similar optimizations that can be applied on the GPU side" is a fascinating question. There are actually two communities trying to write the fastest mandelbrot renderer with very different definitions of fastest. Weirdly, the two communities basically don't acknowledge each other- I don't mean there's animosity, I mean they seem to just not show up in each other's google searches- total lack of mention. On one side is "What is the fastest way to render the mandelbrot set in floating point precision," such as the parent article. This is where you might see tricks like filling in boxes, directly checking for repeating cycles etc, but advanced attempts always seem to get into writing assembly and carefully managing CPU multipliers- The OP's article is a wonderful example of optimizing in this community. I don't really know of many high-man-hour attempts from this branch that use GPU rendering, since it's more about delving into the CPU's details and being faster than other codebases in the same constraints than about achieving a real world task that would otherwise be too slow. The real center of this community is https://benchmarksgame-team.pages.debian.net/benchmarksgame/performance/mandelbrot.html https://benchmarksgame-team.pages.debian.net/benchmarksgame/..., where not only is GPU banned because it would be missing the point, but optimizations like checking whether four corners are all the same color are also verboten, as it's all about how fast the inner loop can spin. The second community is trying to go deep, deep into the mandelbrot set- way past where it makes sense to represent the location of a pixel as a floating point number. Here speed is a practical concern. Images like "Evolution of Trees" https://www.deviantart.com/dinkydauset/art/Evolution-of-trees-402876071 https://www.deviantart.com/dinkydauset/art/Evolution-of-tree... used to take weeks to render, and it's been a no-holds-barred brawl to reduce that time- of course the GPU and alternative algorithms are allowed. The fastest codebases that I know of are FractalShark https://github.com/mattsaccount364/FractalShark https://github.com/mattsaccount364/FractalShark and Kallas Fraktaler 3 https://mathr.co.uk/kf/kf.html#top https://mathr.co.uk/kf/kf.html#top and both have CUDA kernels, but I don't think that these kernels have been optimized to the point of finding and fixing sheep races. The reason is twofold- first, the algorithmic progress has been lightning fast in recent years. The biggest breakthrough, bilinear approximation with rebasing, was discovered in 2022- and it (roughly, the actual runtime is currently unknown) takes the time to calculate an image from O(num pixels * num iterations) to O((num pixels * log (num iterations)) + num iterations). As a result, there just hasn't been time since 2022 to optimize the constant factors to the level in OP, especially now that the algorithm is so much more complicated than multiply, add, repeat. Second, with just lightly optimized CUDA for the num pixels log num iteration term, the overall runtime is currently dominated by the term with no dependency on the number of pixels. No one knows how to parallelize that term (at least not to many cores- fractal shark can use three cores I believe), so improvements in the CUDA kernel don't translate to unlocking deeper renders. The "render deep images fast" community mostly hangs out on fractalforums.org- I highly recommend swinging by, it's a wonderful and active forum (note- heavily moderated to stay on topic- makes hackernews look like the wild west)
alkonaut 3y ago
Cool article. Does anyone have any pointers to what the state of the art is for arbitrary precision? Is that even viable or will it just be a boring extremely slow exercise long before you reach the zoom levels where 64bit float breaks down?
- deleted 3y ago
  [deleted]
- dahart 3y ago
  64 bit doubles break down pretty quickly when zooming into the Mandelbrot: you can zoom in by 2x about 64 times. Doubles only get you twice the zoom level as 32 bit floats, meaning with fp32, you can zoom 2x about 32 times. I feel like using doubles makes zooming a little better but isn’t all that meaningful of an improvement over floats. You can use arbitrary precision, and yes it will be an exercise in patience (both in coding and in waiting for Mandel to render), but there are some really crazy tricks for deep zooming into Mandel that don’t use arbitrary precision, tricks that I know of but don’t understand yet. I’ve seen hints on some of the Mandel deep-zoom videos on YouTube, and also on ShaderToy. BTW I wrote a 128 bit quad-float emulator for ShaderToy, and it seems that it’s somewhere around 100x slower than using floats. It was a very fun exercise.
- QuadmasterXLII 3y ago
  The state of the art for arbitrary precision is in rapid flux right now. The naive approach is O(pixels * iterations * bits_precision^2) which is just hopeless. K. L. Martin[1] took it to O(iterations * bits_precision^2 + pixels * iterations) which was a massive improvement. I explain my webgl implementation of that algorithm here [2] and I think it's currently the fastest web implementation for depths up to 2000 bits. However, in 2022 Zhouran [3] dramatically improved the constant factor for the first term, then made a breakthrough in overall complexity, taking it to O(pixels * log(iterations) + iterations * bits_precision^2) which blows the approach I used out of the water for deep images. It's so recent that there's not really a blog post properly explaining it anywhere, you just have to sort through the forum thread., There is a webassembly implementation at https://fraktaler.mathr.co.uk/live/latest/ https://fraktaler.mathr.co.uk/live/latest/ with public source. 1 http://www.science.eclipse.co.uk/sft_maths.pdf http://www.science.eclipse.co.uk/sft_maths.pdf 2 https://www.hgreer.com/JavascriptMandelbrot/ https://www.hgreer.com/JavascriptMandelbrot/ 3 https://fractalforums.org/fractal-mathematics-and-new-theories/28/another-solution-to-perturbation-glitches/4360/msg29835#msg29835 https://fractalforums.org/fractal-mathematics-and-new-theori...
camel-cdr 3y ago
Somewhat related, ispc also uses the mandelbrot set as a motivating example [0], they get a 5x speedup from writing the code in almost the same way and using ispc. C compilers aren't good at autovectorizing outer loops, which is required here. I recently learend however that intels icx c compiler, based on llvm, can now autovectorize the outer loop [1]. Upstream clang sadly still can't, same goes for gcc. [0] https://ispc.github.io/example.html https://ispc.github.io/example.html [1] https://godbolt.org/z/jbn1robKz https://godbolt.org/z/jbn1robKz
thanatos519 3y ago
There's also the not-a-CPU-optmization where we skip iteration entirely for the main cardioid and circle ... not in C but here's the shader-esque version: ac = vec2(x,y); c2 = dotp(ac,ac); if ((256.0c2c2-96.0c2+32.0x-3.0>0.0)&&(16.0(c2+2.0x+1.0)-1.0>0.0)) { ... do the iterations ... } else { ... it's in the Mandelbrot set ...}
amadio 3y ago
In VecCore (a small C++ SIMD abstraction library on top of Vc and std::simd), I created some simple examples to show how to use the library to optimize code using SIMD in a somewhat generic way. You can find it on GitHub at https://github.com/root-project/veccore https://github.com/root-project/veccore I have examples for Julia sets and the Mandelbrot set, including an implementation with AVX2 intrinsics. These days with std::simd more widely available there's less of a reason to use VecCore, but the examples may still be educational enough. I chose Julia sets and the Mandelbrot since they are perfect examples of simple problems that compilers fail to vectorize on their own.

74 points by RoyWard 3y ago

Could you drop a link on your site to https://www.orange-kiwi.com/index.xml https://www.orange-kiwi.com/index.xml? Hugo generates that RSS feed by default, but you haven't linked it anywhere. Just a friendly heads up :)

tsukurimashou 3y ago
<link rel="alternate" type="application/rss+xml" title="Orange Kiwi" href="https://www.orange-kiwi.com/index.xml" /> yes just add this inside your <head>

deleted 3y ago

[deleted]

dist1ll 3y ago

Cool article. I think what the author calls "sheep race" is known as loop rotation. This allows you to mask off latency, at the cost of wasted computation at the tail of the loop.

RoyWard 3y ago
Do you have a reference for that? I might be looking in different places, but my reading suggests that loop rotation is where a conditional is moved to the end to avoid an unconditional branch, sometimes by the compiler. "Sheep race" with it's wasted computation is a little different - it's not a transformation it would be legal for a compiler to make, and possibly doesn't have enough use cases to really have its own well known name. Oddly enough, I've noticed that the compiler generated code will sometimes not do loop rotation, and and will have a conditional branch in the middle and unconditional at the end. I suspect this is because (a) unconditional branches are often effectively free on modern processors (if all the instructions are cached appropriately), and (b) there are all sorts of subtle instruction layout and ordering that a good compiler does to make things map efficiently to the underlying microarchitecture that means that non-intuitive orderings sometimes have a better cost function. I did find out that what I called "interleaving" is more properly called "software pipelining".

selimthegrim 3y ago

Are the streams problems good concurrency practice?

crazygringo 3y ago

Wow, this article opened my eyes to what is possible in optimizing code for a CPU. I was expecting this article to rehash the usual algorithmic tricks for speeding up Mandelbrot set calculation -- filling in boxes where the four corners already have the same color, checking for repeating iteration cycles to fill in black, etc. But I was very surprised to learn about how important CPU pipelines are, and therefore about "sheep race optimization" and interleaving, and then about the possibility of using vectorization on CPU. And in the end I'm utterly astonished that adding everything up results in the algorithm being calculated eight times faster. I would never have dreamed that such a level of speed improvement would be possible on such a simple, straightforward mathematical algorithm. Of course, in the real world everybody today is going to calculate the Mandelbrot set using the GPU (as the author notes -- this is just an exercise). But now it has me tremendously curious as to whether there are similar optimizations that can be applied on the GPU side? And especially whether those are accessible if done in WebGL/WebGPU? There are quite a number of Mandelbrot set generators available in the browser that use the GPU for calculations (e.g. [1]). [1] http://mandelset.ru/ http://mandelset.ru/ - https://github.com/iskandarov-egor/mandelset https://github.com/iskandarov-egor/mandelset

QuadmasterXLII 3y ago
"Are there similar optimizations that can be applied on the GPU side" is a fascinating question. There are actually two communities trying to write the fastest mandelbrot renderer with very different definitions of fastest. Weirdly, the two communities basically don't acknowledge each other- I don't mean there's animosity, I mean they seem to just not show up in each other's google searches- total lack of mention. On one side is "What is the fastest way to render the mandelbrot set in floating point precision," such as the parent article. This is where you might see tricks like filling in boxes, directly checking for repeating cycles etc, but advanced attempts always seem to get into writing assembly and carefully managing CPU multipliers- The OP's article is a wonderful example of optimizing in this community. I don't really know of many high-man-hour attempts from this branch that use GPU rendering, since it's more about delving into the CPU's details and being faster than other codebases in the same constraints than about achieving a real world task that would otherwise be too slow. The real center of this community is https://benchmarksgame-team.pages.debian.net/benchmarksgame/performance/mandelbrot.html https://benchmarksgame-team.pages.debian.net/benchmarksgame/..., where not only is GPU banned because it would be missing the point, but optimizations like checking whether four corners are all the same color are also verboten, as it's all about how fast the inner loop can spin. The second community is trying to go deep, deep into the mandelbrot set- way past where it makes sense to represent the location of a pixel as a floating point number. Here speed is a practical concern. Images like "Evolution of Trees" https://www.deviantart.com/dinkydauset/art/Evolution-of-trees-402876071 https://www.deviantart.com/dinkydauset/art/Evolution-of-tree... used to take weeks to render, and it's been a no-holds-barred brawl to reduce that time- of course the GPU and alternative algorithms are allowed. The fastest codebases that I know of are FractalShark https://github.com/mattsaccount364/FractalShark https://github.com/mattsaccount364/FractalShark and Kallas Fraktaler 3 https://mathr.co.uk/kf/kf.html#top https://mathr.co.uk/kf/kf.html#top and both have CUDA kernels, but I don't think that these kernels have been optimized to the point of finding and fixing sheep races. The reason is twofold- first, the algorithmic progress has been lightning fast in recent years. The biggest breakthrough, bilinear approximation with rebasing, was discovered in 2022- and it (roughly, the actual runtime is currently unknown) takes the time to calculate an image from O(num pixels * num iterations) to O((num pixels * log (num iterations)) + num iterations). As a result, there just hasn't been time since 2022 to optimize the constant factors to the level in OP, especially now that the algorithm is so much more complicated than multiply, add, repeat. Second, with just lightly optimized CUDA for the num pixels log num iteration term, the overall runtime is currently dominated by the term with no dependency on the number of pixels. No one knows how to parallelize that term (at least not to many cores- fractal shark can use three cores I believe), so improvements in the CUDA kernel don't translate to unlocking deeper renders. The "render deep images fast" community mostly hangs out on fractalforums.org- I highly recommend swinging by, it's a wonderful and active forum (note- heavily moderated to stay on topic- makes hackernews look like the wild west)

alkonaut 3y ago

Cool article. Does anyone have any pointers to what the state of the art is for arbitrary precision? Is that even viable or will it just be a boring extremely slow exercise long before you reach the zoom levels where 64bit float breaks down?

deleted 3y ago
[deleted]
dahart 3y ago
64 bit doubles break down pretty quickly when zooming into the Mandelbrot: you can zoom in by 2x about 64 times. Doubles only get you twice the zoom level as 32 bit floats, meaning with fp32, you can zoom 2x about 32 times. I feel like using doubles makes zooming a little better but isn’t all that meaningful of an improvement over floats. You can use arbitrary precision, and yes it will be an exercise in patience (both in coding and in waiting for Mandel to render), but there are some really crazy tricks for deep zooming into Mandel that don’t use arbitrary precision, tricks that I know of but don’t understand yet. I’ve seen hints on some of the Mandel deep-zoom videos on YouTube, and also on ShaderToy. BTW I wrote a 128 bit quad-float emulator for ShaderToy, and it seems that it’s somewhere around 100x slower than using floats. It was a very fun exercise.
QuadmasterXLII 3y ago
The state of the art for arbitrary precision is in rapid flux right now. The naive approach is O(pixels * iterations * bits_precision^2) which is just hopeless. K. L. Martin[1] took it to O(iterations * bits_precision^2 + pixels * iterations) which was a massive improvement. I explain my webgl implementation of that algorithm here [2] and I think it's currently the fastest web implementation for depths up to 2000 bits. However, in 2022 Zhouran [3] dramatically improved the constant factor for the first term, then made a breakthrough in overall complexity, taking it to O(pixels * log(iterations) + iterations * bits_precision^2) which blows the approach I used out of the water for deep images. It's so recent that there's not really a blog post properly explaining it anywhere, you just have to sort through the forum thread., There is a webassembly implementation at https://fraktaler.mathr.co.uk/live/latest/ https://fraktaler.mathr.co.uk/live/latest/ with public source. 1 http://www.science.eclipse.co.uk/sft_maths.pdf http://www.science.eclipse.co.uk/sft_maths.pdf 2 https://www.hgreer.com/JavascriptMandelbrot/ https://www.hgreer.com/JavascriptMandelbrot/ 3 https://fractalforums.org/fractal-mathematics-and-new-theories/28/another-solution-to-perturbation-glitches/4360/msg29835#msg29835 https://fractalforums.org/fractal-mathematics-and-new-theori...

camel-cdr 3y ago

Somewhat related, ispc also uses the mandelbrot set as a motivating example [0], they get a 5x speedup from writing the code in almost the same way and using ispc. C compilers aren't good at autovectorizing outer loops, which is required here. I recently learend however that intels icx c compiler, based on llvm, can now autovectorize the outer loop [1]. Upstream clang sadly still can't, same goes for gcc. [0] https://ispc.github.io/example.html https://ispc.github.io/example.html [1] https://godbolt.org/z/jbn1robKz https://godbolt.org/z/jbn1robKz

thanatos519 3y ago

There's also the not-a-CPU-optmization where we skip iteration entirely for the main cardioid and circle ... not in C but here's the shader-esque version: ac = vec2(x,y); c2 = dotp(ac,ac); if ((256.0c2c2-96.0c2+32.0x-3.0>0.0)&&(16.0(c2+2.0x+1.0)-1.0>0.0)) { ... do the iterations ... } else { ... it's in the Mandelbrot set ...}

amadio 3y ago

In VecCore (a small C++ SIMD abstraction library on top of Vc and std::simd), I created some simple examples to show how to use the library to optimize code using SIMD in a somewhat generic way. You can find it on GitHub at https://github.com/root-project/veccore https://github.com/root-project/veccore I have examples for Julia sets and the Mandelbrot set, including an implementation with AVX2 intrinsics. These days with std::simd more widely available there's less of a reason to use VecCore, but the examples may still be educational enough. I chose Julia sets and the Mandelbrot since they are perfect examples of simple problems that compilers fail to vectorize on their own.