Loop unrolling is complexicated. Usual knowledge is that the compiler has way more knowledge than you might have about source code and the target machine, and so you should just write a loop, and let the compiler decide what to actually do with it. This is a good rule of thumb, most of the time.

The other day I made one GPU compute shader in Blender 10 times (ten times!) faster by unrolling a loop though. (ten times on one GPU… on another GPU it is two times faster, but still!)

Kuwahara image filter

Blender has a Kuwahara compositor node, that produces a “painterly” look from an input image. Sometimes it is used as a building block for more complex painterly effects.

Fun fact! Original Kuwahara algorithm was developed for biomedical image processing needs, to smooth out noise while preserving boundaries. Kuwahara, Hachimura, Eiho, Kinoshita “Processing of RI-Angiocardiographic Images”, 1976. Much later, Kyprianidis, Kang, Döllner “Image and Video Abstraction by Anisotropic Kuwahara Filtering” (2009) extended it to look better and the goal was non-photorealistic image & video processing.

Anisotropic Kuwahara filter is quite expensive to calculate, and most of the cost comes from the fact that each and every output pixel must look at an oriented ellipse around itself (ellipse orientation & shape changes from pixel to pixel), sample input pixels inside of it, assign the inputs to eight sectors, calculate sector statistics, and then blend their average colors, preferring sectors with lower variance. The sectors are also actually overlapping with each other, ugh.

But, it produces quite nice “painterly” results, here’s a result that we can consider to be ground truth because the dog’s name is True:

The problem

Anyway, I was looking at it for some small optimizations. On my PC (Ryzen 5950X / RTX 3080Ti), doing Kuwahara with size 8 on a 2048x1536 input image took 408ms on the CPU, and 19ms on the GPU. So the GPU is about 20x faster at this, which sounds fine.

But then! On an Apple M4 Max device, the CPU took 385ms (faster than this Ryzen, eh), yet the GPU took a whopping 220ms. Like, what? One could go “well obviously this Apple GPU is slower than that nvidia GPU! And GPU path on the Apple device is still faster than CPU, so good” and move on. But the ratio felt way off; why the first machine would have 20x performance advantage by using GPU, whereas the second one would not even be twice as fast? And also, the M4 Max GPU might be slower than RTX 3080Ti, but it should not be ten times slower.

Let’s investigate a Metal shader

Apple Xcode can do a frame capture of an application that uses Metal graphics API, and then you can use it to do debugging, performance insights and so on. I’m using Xcode 26.3 on macOS 15.7.3. Capture a frame with performance analysis enabled, and click Show Performance button:

The Top Shaders view shows the most expensive shaders, where we find our main Kuwahara compute shader. It uses 192 registers (a lot!), and even spills 64 bytes since it’s kinda out of registers (that is not great at all):

Clicking on Counters tab makes the whole computer almost unresponsive for several minutes, while Xcode is really busy doing something. Then it shows that our compute shader only has 28% occupancy (not great), but more curiously, the instruction mix is like:

• ALU float instructions: 12% (huh? that low? source code was full of math!)
• ALU half instructions: 0% (yes, Blender today does not use FP16 math in shaders… perhaps it should, but that’s for another day)
• Conditional, Integer and Boolean instructions: 87%. What?!

With our compute shader still selected in Overview tab, right click on it and pick Reveal in Cost Graph, and it shows a flamegraph of the shader code & function calls (pretty much everything in just main entry point), but more importantly, it shows the shader source code with approximate per-line costs below. Scroll to the expensive part:

This is the inner shader loop, where for each image sample we check which sector(s) it is in, weigh them, calculate their statistics and so on. It would make sense that this is the heaviest part of the shader, since it is the innermost loop after all. But hovering on each of the line-attribution-bubbles shows something like this:

Wait what? Why would a line like weighted_mean_of_color_of_sectors[upper_index] += upper_color * weight be mostly “select” instructions, and not just math?!

At least this version of Xcode just stops there. It says there are a ton of “select” instructions, but does not actually show them. What is going on there? Why? No idea, but it would be really nice if Apple had tools that could enable us to drill deeper. Since Xcode can not help us further, we’ll have to go heavy metal 🤘

Dougall Johnson has a project with Apple M1 (“G13”) GPU reverse engineered, with extensive documentation and disassembler and other tools at github.com/dougallj/applegpu. It is not exactly the same GPU that I am on (which is M4), but this will have to do.

After some goofing around, something like this seems to achieve a “GPU disassembly for M1 GPU” from my Metal source code:

1. Compile Metal source code into a metal library: xcrun -sdk macosx metal metal_source.metal -o test.metallib
2. Create a test.mtlp-json file (this describes the pipeline state) with contents like:

   {"pipelines":{"compute_pipelines":[{
   "compute_function":"_compositor_kuwahara_anisotropic_constant_size_comp",
   "threadgroup_size_is_multiple_of_thread_execution_width":true,
   "max_total_threads_per_threadgroup":1024}]}}    
   

   and then compile to M1 GPU (applegpu_g13g) code with xcrun metal-tt -arch applegpu_g13g -o test_g13g.gpubin test.mtlp-json test.metallib
3. This produces a .gpubin file which is a Mach-O file, and we are interested in the __TEXT,__compute section. That section itself is another Mach-O file! And from that file, we need the __TEXT,__text section. With that extracted into a separate file like g13g.bin, we could finally disassemble it.
4. But the disassemble.py tool stops on the first stop instruction, which is placed in some sort of shader preamble. Disable this stopping to get full shader disassembled:

   cd applegpu
   python3 -c "
   import disassemble as d
   d.STOP_ON_STOP = False
   d.disassemble(open('../g13g.bin','rb').read())" > ../g13g.asm
   

And yes, majority of instructions in the shader disassembly are icmpsel. There’s also a bunch of stack_store and stack_load, which seem like that’s the “spilled bytes” that Xcode was pointing out. Here’s a snippet of actual instruction sequence from the GPU shader disassembly:

mov              r9, r51.cache
icmpsel          seq, r51, r1l.cache, 2, r7.cache, r34.discard
icmpsel          seq, r18, r1l, 3, r7, r18
stack_store      i16, 1, 0, xy, 4, r18l_r18h, 88, 0
icmpsel          seq, r16, r1l, 4, r7, r44.discard
stack_store      i16, 1, 0, xy, 4, r16l_r16h, 84, 0
icmpsel          seq, r7.cache, r0h.cache, 1, r14.cache, r12.cache
icmpsel          seq, r14.cache, r0h.cache, 1, r39.discard, r6
icmpsel          seq, r7.cache, r0h.cache, 2, r51, r7.cache
icmpsel          seq, r14.cache, r0h.cache, 2, r24.discard, r14.cache
icmpsel          seq, r7.cache, r0h.cache, 3, r18, r7.cache
icmpsel          seq, r14.cache, r0h.cache, 3, r13, r14.cache
icmpsel          seq, r7.cache, r0h.cache, 4, r16, r7.cache
icmpsel          seq, r14.cache, r0h.cache, 4, r5, r14.cache
icmpsel          seq, r7.cache, r0h.cache, 5, r38, r7.cache
icmpsel          seq, r14.cache, r0h.cache, 5, r3.cache, r14.cache
icmpsel          seq, r7.cache, r0h.cache, 6, r20, r7.cache
icmpsel          seq, r14.cache, r0h.cache, 6, r11, r14.cache
icmpsel          seq, r7.cache, r0h.cache, 7, r19, r7.cache
icmpsel          seq, r14.cache, r0h.cache, 7, r8, r14.cache
fadd32           r7.cache, r7.cache, r31
fmadd32          r50, r32, r29, r14 

Out of 22 instructions, only two do some sort of math! The rest are integer comparison-selects and some stack data movement. What seems to be going on, is that the final GPU shader compiler did not unroll this loop (every loop iteration has the array index [i] that is just a variable), and shader variables like float sector_weights[8], float4 weighted_mean_of_color_of_sectors[8] and float4 weighted_mean_of_squared_color_of_sectors[8] become actual arrays inside the GPU registers. At least this GPU cannot index registers! So in order to do any array index the compiler emits a sequence like if i==0 use this, else if i==1 use that, else if i==2 use this, .... All of that expands to about 200 integer select instructions per loop iteration.

Madness!

Let’s unroll the loop

Blender’s shaders are written in a curious shading language tentatively called “BSL” (see BCON26 talk about it, “How Not to Build a Shading Language”). For explicit loop unrolling it has a C++-style attribute [[unroll]], which literally copy-pastes the loop body before passing it to a platform compiler.

for (int k = 0; k < 8 /* number_of_sectors */; k++) [[unroll]] {
  float weight = sector_weights[k] * radial_gaussian_weight;
  sum_of_weights_of_sectors[k % (number_of_sectors / 2)] += weight;

  int upper_index = k;
  weighted_mean_of_color_of_sectors[upper_index] += upper_color * weight;
  weighted_mean_of_squared_color_of_sectors[upper_index] += upper_color_squared * weight;

  int lower_index = (k + number_of_sectors / 2) % number_of_sectors;
  weighted_mean_of_color_of_sectors[lower_index] += lower_color * weight;
  weighted_mean_of_squared_color_of_sectors[lower_index] += lower_color_squared * weight;
}

This makes the shader run 10 TIMES faster.

It was taking 220ms before, and now it runs in 20.4ms. The shader instruction mix went from (12% float math, 87% conditional and integer) to (79% float math, 17% conditional and integer), which feels way more sensible. It still uses 192 registers and spills (64 bytes -> 48 bytes), but that’s a topic for another day. And the (unrolled) inner loops are now all just math instructions:

Okay, but what about other GPUs?

Fair question! I only have two GPUs around here, the other one being RTX 3080Ti. Which, funnily enough, is another GPU where the GPU vendor does not easily allow you to see “what is the actual GPU assembly” bit :)

Nsight Graphics allows me to debug and profile graphics applications, but for example when using Vulkan, it can show me which of the SPIR-V assembly instructions were costly, but how they map into the underlying GPU instructions, is anyone’s guess.

In this particular case, speedrun of Kuwahara through Nsight (2026.3 version): Start Activity, pick “GPU Trace Profiler” option, turn on “Real-Time Shader Profiler”, launch the application, do a capture. Find the expensive compute shader dispatch in the timeline (it takes 17.8ms here):

In lower right corner, under Shader Pipelines section, it displays some data about this compute pipeline:

It can run 8 warps, uses 204 registers (198 live registers), average warp latency is 134K, instruction mix: 44% FMA FP32, 36% data movement.

With unrolled loop, now our compute dispatch takes only 9.5ms:

And it can run 16 warps, uses 113 registers (103 live registers), average warp latency is 45K, instruction mix: 60% FMA FP32, 13% data movement.

So yes even on RTX 3080Ti this particular shader gets twice as fast by just unrolling the inner loop. Half the registers results in twice the occupancy here, but TBH I did not dig much deeper into why exactly the register usage has decreased.

What’s the lesson?

Yes the loop unrolling should be left to the compiler… except when it should not :)

That by itself is not very useful, so perhaps a better lesson is - it helps to have at least a vague understanding of what sort of performance should be expected. This way you can go “wait a minute, this feels way too slow” and investigate. In this case, I wanted to fix seemingly very bad compute shader performance on an Apple GPU (and sped it up ten times), but a side effect was that on an NVIDIA GPU the shader also got twice as fast.

It helps if the performance analysis tools can point out some sort of data that makes you go “huh wait, this should not be here”, like in this case, majority of shader instructions being integer compare/selects. It would be even better if these tools could actually show more data; for Apple GPUs in Xcode frame capture and NVIDIA GPUs in Nsight, I’d like to see GPU assembly. One can dream!

And now I’ll go and ship these Kuwahara node performance improvements in Blender (PR 164617, just merged for upcoming Blender 5.3). There’s probably more work to do in that area, but that’s for later.

Almost a year ago, while testing some Blender compositor related things, I noticed that their “Fast Blur” mode is not very fast in some cases, and I started to wonder maybe Blender should get a Dual Kawase blur mode, but that had limitations, and after a whole bunch of poofing around I assembled a bunch of existing blurry ideas into Smol Gaussian.

So! As everyone knows, doing a Gaussian blur in the most naïve way is a cost that grows up as square of the blur radius. But luckily enough, Gaussian is separable, so you can do horizontal & vertical blurs separately, and now the cost scales linearly with radius. That is perfectly fine for blurring just a little bit, but if you need really large blurs, that becomes slow.

There are many ways to make blur scale to larger radii, FrostKiwi has a really good blog post about them: Video Game Blurs (and how the best one works), and the “best one” that the post ends with is Dual Kawase.

Dual Kawase is really simple to implement, runs really fast, and is very nice! Except if you 1) need a smoothly varying blur radius, or 2) need different amounts of horizontal & vertical blur. In my case, I needed both. So I kept on tinkering with trying to extend Dual Kawase, and while I got something that kinda works, I found another way to blur that is much closer to Gaussian shape, and feels nicer when blur radius is animated. Or to phrase it differently, I have discovered what a ton of other people have done as well, which is “downsample, blur, upsample”. There’s quite some details to get through though :)

If your daily reading limit is already up, you can stop here and go check an interactive WebGPU blur playground. Otherwise, continue reading!

In the videos below, blur radius is animated (increase X&Y together from 5 to 1000, decrease X&Y separately), and the input image is rotated during the animation loop. The input image has very bright objects (linear intensities 500 and 50), here’s the input with additional glare to show “hey I’m bright” (actual image used for tests has no glare).

Smol Gaussian

Here’s my Smol Gaussian: downsample to a working resolution (which could be different per axis), apply a small separable Gaussian there, reconstruct to full resolution image.

Here’s how it works in imaginary pseudocode:

Image smol_gaussian(Image input, float2 radius)
{
    float2 sigma = max(radius, 0) / 3;
    int2 orig_size = input.size;
    int2 working_size = orig_size;
    Image image = input;

    // Downsample to the working resolution.
    while (true)
    {
        float2 scale = orig_size / float2(working_size);
        bool2 reduce = (sigma / scale >= 6) && (working_size > 1);
        if (!any(reduce))
            break;
        int2 next_size = select(working_size, ceil(working_size / 2.0), reduce);

        // Integrate source pixel areas, taking care of odd dimensions.
        // Keep one pixel border on reduced axes to preserve input edges.
        // Two consecutive exact halves can be merged into single 4x reduction.
        image = downsample_with_border_preserve(image, next_size);
        working_size = next_size;
    }

    for (axis in {X, Y})
    {
        float scale = float(orig_size[axis]) / working_size[axis];
        float down_variance = (scale * scale - 1) / 12;
        float up_variance = scale == 1 ? 0 : scale * scale *
            (scale > 2 ? 1.0/3 : scale == 2 ? 3.0/16 : 1.0/6);
        // Variance estimates are in original pixels; convert to working pixels.
        float residual_sigma = sqrt(max(0,
            sigma[axis] * sigma[axis] - down_variance - up_variance)) / scale;

        // Regular normalized 1D gaussian pass, extending to 4*sigma, smoothly
        // tapered to zero between 3*sigma and 4*sigma. Pair adjacent weights
        // into bilinear texture samples; skip axes with negligible sigma.
        image = gaussian_pass(image, axis, residual_sigma);
    }

    // Per axis: bilinear up to 2x enlargement, cubic B-spline beyond 2x.
    return reconstruct(image, orig_size);
}

The algorithm is similar to Skia’s GPU Gaussian blur as of 2026 Sep (FilterResult::Builder::blur and FilterResult::rescale in SkImageFilterTypes.cpp) - independent X/Y scaling, texture samples placed to use bilinear filtering, one pixel border on downsampled images that preserve the original image edges (this way bright interiors do not overbright the result at large radius).

Additions compared to Skia blur are:

• When downsampling odd-sized image, we do more correct pixel area integration, so that isolated bright pixels do not flicker with sampling phase shift.
• Downsampled levels are ceil-halved sizes, and downsampling kicks in at sigma=6 (so in practice, blur sizes 18, 36, 72, 144, … switch to new level). Skia instead scales continuously, to keep working sigma under 4.
• When doing the final Gaussian blur, we take into account the blur introduced by downsampling and later reconstruction, i.e. subtract their variance from the blur kernel.
• The final Gaussian kernel extends to sigma=4 (i.e. not truncated at sigma=3), and weights in sigma 3..4 region are tapered to reach zero. This helps to reduce the “blur is cut off” with very bright highlights, and looks better when radius is animated and number of taps changes.
• Final reconstruction to full image size uses a cubic B-spline (fully positive kernel, so no ringing) instead of bilinear, beyond 2x enlargement. This helps to avoid slope discontinuities of bilinear.
• Pairs of exact 2x reductions are done as single 4x reduction as an optimization.

Related work for all of this:

• Skia SkImageFilterTypes.cpp, as mentioned above.
• Fabian Giesen, “Gaussian blur kernels” on gdalgorithms list (2009): low-pass filter, blur at lower resolution, then upsample. He suggests bilinear upsampling is good enough, however in my tests very bright (HDR) blurred objects show bilinear “slope steps”, so I used cubic.
• Intel, “An Investigation of Fast Real-Time GPU-Based Image Blur Algorithms” (2014) has “Working in Lower Resolution” section, but with not much details at when you would switch to it, or how to handle animated blur radius.
• GPU Gems 2, Chapter 20, “Fast Third-Order Texture Filtering” has a trick for cubic B-spline filtering using bilinear samples, which is what I used here in reconstruction part.

All of this is not particularly novel or insightful, I consider it to be perhaps a handful of small quality tweaks compared to Skia blur.

A discarded idea: crossfading working resolutions.

I tried blending neighboring resolutions to hide level switches. Idea was this: when the resolution that does the final blur changes, in theory you could have a visible “jump” if blur radius is animated.

So if resolution switch happens at sigma=6, then starting at sigma=5 already, do both the current resolution and the next resolution, and blend between them using a smoothstep curve. Both evaluations target the same final blur amount, just at different grid, and so each of them uses different amount of gaussian taps.

This however cost performance, since now during transition regions (blur sizes 15..18, 30..36, 70..72 etc.) there are two Gaussian blurs performed and their result is blended. If the X/Y radii are different, we might need to evaluate three Gaussian blurs in fact, and blend between them.

In my testing, this brought pretty much no visual difference, but cost in performance and made the implementation more complex. So eventually this was discarded.

Dual Kawase

The playground also implements an extended version of Dual Kawase. See Marius Bjørge, Bandwidth-Efficient Rendering (SIGGRAPH 2015). The original uses equal horizontal/vertical blur amounts, and only supports a discrete “number of blur pyramid levels” control, not a continuous “blur radius” setting.

For arbitrary blur sizes: blend between neighboring discrete blur levels, similar to obs-composite-blur. Blend using the fractional position between steps remapped with t * (2 + t) / 3, this makes it feel a bit nicer than just a linear blend.

For independent X/Y blur radii: when the smaller blur radius is reached, we stop further reductions along that axis. For the between-levels blend above, we might need to blend between three different blurred results.

This works, but does not “feel great” to me. Smoothly animated blur radius does not “feel” smooth on HDR highlights, due to how blending happens between discrete blur levels.

Skia Gaussian

This is what conceptually is closest to “Smol Gaussian”. Implementation in my playground is just a WebGPU re-implementation of relevant parts of Skia SkImageFilterTypes.cpp, as it was in revision 15a9437eec87 (2026 Sep). Basic algorithm is:

• Each axis independently downscales to a working sigma <= 4. Intermediate steps halve the scale, and the last step uses the remaining fractional scale.
• A one pixel border around intermediate steps preserves clamped edge colors.
• Final Gaussian pass has radius of ceil(sigma * 3). A single direct convolution, or separable bilinear-paired filter samples is used depending on size.
• Final result is bilinearly upscaled to original resolution.

Skia blur feels just fine on regular LDR content, but on very bright HDR highlights, the aliasing and wobbling are quite apparent.

“Fast Gaussian” (from Blender 5.2)

Blender’s compositor blur node got a “Fast Gaussian” mode back in 2008 (v2.46) in 2a2453d3, which however built upon earlier implemented IIR_Gauss functionality for defocus blur node (2006, v2.43, commit e61dec07).

This builds upon “Recursive Gaussian Filtering”, which if you’re just a programmer but not familiar with signal processing terminology, is quite a confusing name. You’d think a recursive gaussian would be something about building like several smaller versions and somehow combining them, right? Haha nope, not at all, “recursive filter” in signal processing just means that the filter uses some of the previous outputs.

Anyway, the “fast” part is due to filter construction that is basically the same cost, no matter the blur radius. Original code in Blender seemingly was built on one such algorithm, from Young, van Vliet & van Ginkel, “Recursive Gabor Filtering” (2000) paper. Many years later, in 2024 (blender 4.2.0, commit 382131fe) Omar remade this algorithm to have both CPU and GPU code paths, and to avoid double precision. And it was built on a handful of papers, curiously enough an earlier paper by Young, van Vliet et al. “Recursive Gaussian derivative filters” (1998), another paper Deriche “Recursively implementing the Gaussian and its derivatives” (1993), and some more.

Anyhoo, in this project there’s a WebGPU re-implementation of Blender 5.2 “Fast Gaussian” state (mostly recursive_gaussian_blur.cc), which is fourth order Deriche formulation for radius under 96, and Van Vliet formulation for larger radius.

The algorithms are more or less constant work independent of the blur radius, which is very nice. However, they are also from 25+ years ago, and are not “embarrassingly parallel” that would fit a GPU (or even a many-core CPU) very well. So despite the name, they may or might not be very fast :)

They also have some ringing artifacts, which are not that much noticeable in regular colors, but with very bright HDR highlights, blurred result can have halos or negative colors, which is not great. See:

“Ryg Blur”

I’m calling this “Ryg Blur” due to Fabian Giesen’s Fast blurs 1 and Fast blurs 2 blog posts (2012), which nicely describe the whole idea, including how exactly to handle fractional samples at the ends, and how trivially that extends to a compute shader implementation.

However the idea itself is “repeated box convolution”, and has been around for ages, e.g. Heckbert “Fun With Gaussians” (1985) talk about it in pages 11-12.

“convolution of a 500x500 image with a 35x35 kernel would take over an hour with 2-D convolution, but only 4 minutes with two 1-D convolutions” where he’s talking about a separable Gaussian filter. Only four minutes to blur a 500x500 image, imagine that! Computers have gotten quite a bit faster, eh.

Anyway, this implementation is basically fixed cost independent of blur size, and uses two bilinear samples for fractional box endpoint updates. Number of iterations control is for how many times this box convolution should be done (1: box filter, 2: tent filter, 3 and up: approaching Gaussian). It is very simple to implement, however not the fastest. Amount of parallelism (parallel over rows or columns) is nowhere near enough to feed modern GPUs, and multiple passes over the full size image incur a lot of memory traffic.

Regular Gaussian blur

There’s also of course a simple separable Gaussian blur, mostly included as a reference. This is one algorithm that becomes impractical at very large blur radii, since the cost scales linearly with radius. This particular implementation cuts off the kernel at sigma=3 (matches behavior of Blender), which is fine for regular image content, but for very bright HDR highlights it makes the blur feel like it “stops” abruptly. Smol Gaussian feels better in this regard!

Performance

All of the below are in Chrome browser. The timings finish any GPU work, start measuring, do the blur, finish GPU work again, end measuring. So this kinda measures “latency” of doing both CPU & GPU work parts.

Apple M4 Max, macOS:

Ryzen 5950X, RTX 3080Ti, Windows:

Core i7-1185G7, Intel Iris Xe, Windows:

Summary:

• Smol Gaussian is about same performance as Skia blur (same on RTX 3080Ti, a bit slower on M4 Max, a bit faster on Intel Xe). Both of them actually get cheaper as blur radius increases! 🤯 This is because the costly part (actual Gaussian convolution) happens on smaller & smaller image, with the rest being quite cheap.
• Dual Kawase is 2x-3x slower than Smol Gaussian, at least in this implementation that blends between levels and tries to handle non-isotropic blur.

More on the blur playground

Again, there’s an interactive blur playground and source code for all of it. Requires WebGPU support (including float32-filterable and float32-blendable).

• Pick any of the six blur modes! For free!
• Drag & drop or browse for images (PNG, JPG and a subset of EXR). test_a_small and test_b_1080p links load two example EXRs.
• Adjustable blur radius, with ability to lock X/Y.
• Can inspect various intermediate textures produced by a blurring algorithm.
• Animate checkbox animates the radius and the source image rotation.
• Render Video button exports animated blur result into a video file, using the current blur mode. The video is encoded using WebCodecs browser functionality, and resized to max 960px. Blur itself is performed at full image resolution.
• Benchmark button tests increasing blur radii with all the selected blur algorithms, and produces a SVG file with the graphs. The result is displayed at the bottom of the page, and can be downloaded too.

Disclosure: I had lollum clankers write javascript & webgpu parts, since I know nothing about this technology. It was both an impressive and quite annoying experience!

Should I try to get this into Blender?

All of this started almost a year ago when I got curious why “Fast Gaussian” does not feel fast in some cases. Well, now I have a WIP pull request with the findings for Blender compositor node (PR 164496). We have not decided what to do with it yet; e.g. should it just replace the existing “Fast Gaussian”, or should it be a new blur mode, or what. That PR also has a CPU implementation (Blender’s compositor today needs to support both GPU and CPU execution), and while some details are different (no paired bilinear sampling on the CPU), it is still several times faster than the current Fast Gaussian, and with no ringing artifacts either.

So, will see how this goes!

Also, if I managed to get something completely wrong about how blurring should work, or there are obvious improvements to do, let me know via email or on mastodon. See ya!

Here’s a bunch of random things that happened in Blender Video Sequence Editor (VSE) since the previous blog post.

For the full picture check out VSE related release notes of Blender 5.1, 5.2 and upcoming 5.3. Below I’m mentioning only the things that I have done, or was somewhat adjacent to the things being done.

Compositor ❤️ VSE

Blender 5.0 added compositor based strip modifiers: in addition to a handful of existing modifiers for visual strips (mostly various color adjustments), you can now make completely custom modifiers that are driven by a compositor node tree. So now you can do anything that Blender compositor nodes allow you to do. Nice!

However, that was only the initial glimpse of the things to come, and it required quite some follow-up work in order for it to be actually useful :)

Compositor based transitions

Strip modifiers work on a single strip, and they alter the image produced by it. “Effects”, on the other hand, can have multiple strips as inputs, and produce an image based on the inputs and often some sort of “time factor”. The most common case is transition effects: two inputs, and a 0..1 range number indicating how far into the transition we are. Blender had only very simple built-in transition strips: “cross fade” and “wipe”.

Now, starting with 5.2, you can have compositor based effects (PR #150694), and one could implement a ton of custom transitions.

Transitions are perhaps the most common, but you can also have effect strips that do not have any inputs, and just generate an image completely procedurally. This would be useful for any sorts of noise, patterns or other procedural textures.

Custom inputs for compositor based modifiers & effects

While having a node tree to drive a modifier or effect is nice and very flexible, it is of limited use if you need to hardcode every tweakable parameter directly in the node tree. It would be way more useful if you could prepare a node tree once, turn it into an asset, and expose whatever parameters a user might want to tweak as custom inputs. Right? Right.

So that’s what Falk implemented for compositor modifiers in Blender 5.2 (PR #156706), and I did the same for compositor effects in upcoming 5.3 (PR #161026).

More performance around compositor

The initial implementation of compositor modifiers in 5.0 was always using the CPU based compositor. This works, but is not the fastest way of doing large amounts of pixel operations :) And so I spent some time improving various aspects of it:

• The obvious one is, let it use the GPU compositor, when the user has selected the GPU compositor (the GPU compositor has been recently made the default option, by the way). PR #154629
• While at it, also allow using half precision data across the GPU compositor. PR #154914
• Do all colorspace conversions “lazily” (right before they are needed) across the whole VSE rendering stack. Previously, some cases were doing too many “back and forth” conversions, especially when scene linear & display colorspaces are intermixed. This helps compositor use cases (the compositor in Blender always uses scene linear space), but also cases like having sequences of EXR images inside VSE, and so on. PR #156281
• While at it, do compositor colorspace conversions on the GPU, when using the GPU compositor. PR #156640
• Automatically detect when the resulting image of a compositor modifier or strip is fully opaque. The initial implementation used to be conservative, and assumed the result might contain transparency. Which also means that visual strips that are fully beneath the current strip could not be “culled away”, since we did not know if the strip is fully opaque or not. Not anymore! Now it detects any non-opaque-alpha presence (via parallel reduction) and properly marks the result as opaque or not. PR #158600
• Did various optimizations for Mask strips. PR #155853, PR #156518, PR #156622.

There’s still more things to do on the performance front. The most obvious ones being: 1) more of the whole VSE rendering stack being able to use the GPU, and 2) hardware accelerated video encoding/decoding. These (and more) are on the 2026 roadmap, but given how fast time seems to be going, we’ll see how much will actually get implemented this year.

More thumbnails for more things

Image, Movie and Sound strips already had thumbnails, but some other strip types did not. One of the unique features of Blender VSE is that you can directly use other scenes as strips in your timeline, without first rendering them to a video or image sequence. This seems to be used all the time at Blender Studio for storyboarding, for example. Alas, scene strips did not have thumbnails. Well, they will in Blender 5.3 (PR #160322)!

Text, Mask and MovieClip strips also learned how to have thumbnails (PR #151042).

Reduce frame drops around video cuts

It is very common to have an input video file be chopped up into smaller pieces while editing:

However, previously each and every strip maintained its own movie decoder object (internally that contains an ffmpeg decoding context, one or more decoded frame buffers, and so on). Now, a single decoder object is fairly heavy, and can consume like a hundred megabytes of memory (depending on video codec and resolution), so you would not want all the video strips in the timeline to have them open at once. However, initializing and closing these movie decoders is also slow (ffmpeg does whatever it needs to do to open & recognize a file, plus there’s a handful of large memory allocations). What the previous code was doing, is that video strips that are playing right now would keep their movie decoders open, but as soon as the playhead moved past them, they would close now-stale decoders.

Kinda makes sense… except for video cuts. Two strips that would be neighboring, and read the same input video file, still had their own completely separate movie decoders. The one that just finished playing would close the decoder, and the one that just started playing would open its own decoder. Which would do all these heavy memory allocations and ffmpeg stream/codec detection work all over again.

Ugh. This meant that (if you had VSE frame prefetching turned off), for interactive playback you would often get frame drops right around video cuts; indicated by gaps in the blue line here:

Now of course, frame prefetching (which is on by default) somewhat hides this. But still, doing this heavy work around each and every cut is wasteful.

Well, not anymore! For Blender 5.3 (PR #162451) I have implemented a shared movie decoder pool/cache. No longer does each strip maintain its own movie decoder; instead they are served from a shared cache, and the cache tries to select the most suitable currently-unused decoder, taking into account how far from the needed frame the last decoded video position is, and so on. This not only fixes the interactive playback frame drops, but also saves quite a bit of final render time, when you have many video cuts.

There have been several attempts in the past at fixing this exact issue; I’m quite happy we finally landed one.

Misc

As part of Code Quality Project 2025 I did a whole bunch of cleanups, refactors and ancient code removals across the VSE codebase. Removing code is the best of those, of course :)

I did some speedups and fixes to proxy building for images (PR #152891) and movies (PR #152949), but nothing to write about. Except ha! I just did! My blog, my rules.

As part of an investigation and fixes to some audio related memory leaks, I completely rewrote the external library (Audaspace) integration into Blender code (PR #153386). I wish this had been about some thread safety fixes, so I could quote a song like “I’m gonna send him to Audaspace, find another race”. Alas, this was about memory leaks, so I can’t.

Oh! And the most important thing. Last time I said I became the “lead” of the whole VSE effort. Well, good news, that only lasted like half a year, and I passed the title & responsibility to John Kiril Swenson, who is doing an excellent job. At some point I should learn that each and every time I try to become a “lead” of anything, that only leads to misery for everyone involved.

okay bye!

A Blender issue #152467 (“File Browser thumbnail cache broken with large amount of images”) reminded me to write this up. This particular issue is a (documented) surprise that when you have a parallel_for in TBB, some of the loop iterations might not execute at all, if the task group gets cancelled.

Similar to C++ exceptions, the effect is “global” - something might throw an exception, and now a completely unrelated part of your code needs to be aware of that possibility, even if you don’t want to. Same with task group cancellation – you might write a parallel_for, and assume that all the loop iterations will execute. That’s what all the code within Blender does, I think :) But! Because some caller of your code way up above might do a task_group.cancel(), now your code needs to be prepared to handle that possibility.

Anyway, all that reminded me of another Blender bug that I was involved with some months ago, which is more curious.

The bug

The reported bug was #143662:

• Blender crashes, while you have a file browser dialog open with “sufficient amount” of thumbnails,
• But only if some thumbnails were freshly generated (i.e. have not been cached previously),
• And only if you had rendered anything with the path-tracing Cycles renderer,
• And only if you have “Persistent Data” Cycles option on,
• The crash does not happen with Address Sanitizer being on.

So that’s… fun.

Part of the possible cause was my change that added more multi-threading to parts of image processing code, some of which gets executed during thumbnail generation. Which means I had to investigate.

What is wrong with this code?

Suppose you have a path-traced renderer (like Blender’s Cycles) that as part of scene initialization does something like this:

parallel_for_each(all_scene_geometries, build_geometry_bvh); // 1.

and each Geometry object contains something like:

struct Geometry {
    BVH bvh; // bounding volume hierarchy object backed by Embree library // 2.
    // ...
};

So far so good. This builds bounding volume hierarchies for all geometries, in parallel. The parallel_for_each is implemented by TBB library, and the BVH data is backed by Embree library. Both well known, production & battle-tested libraries.

Now, in a completely unrelated part of the application, like in a file dialog UI code, you have an on-demand file thumbnail generation code. Thumbnails are cached to disk, but if some of them are not cached, they are rebuilt in the background and saved. While scrolling the dialog with potentially many thumbnails, some of the queued requests might get no longer needed if the visible portion changes drastically. The exact logic is somewhat complex, but essentially it has a queue of “thumbnail generation tasks”, and sometimes decides to cancel a whole group of pending tasks.

So there’s code like this somewhere:

if (some_condition) {
    tbb_thumb_task_group->cancel(); // TBB cancel functionality // 3.
    // ...
}

and somewhere within each thumbnail generation job, there is potential image buffer colorspace conversion or scaling code that has like a:

parallel_for(all_thumbnail_pixels, process_the_pixel); // 4.

Does all of that code make sense? You would think so! And yet it crashes, from innards of oneTBB, when doing the cancel() call. But only sometimes. And only when Address Sanitizer is off.

And only if all of the 1, 2, 3, and 4 points above are present:

• Remove the parallel build of Cycles geometry BVH, i.e. build them sequentially? All good.
• Switch the Cycles BVH to something not backed by Embree? All good.
• Switch the Cycles to not persist scene data across frames/renders? All good.
• Stop doing cancel on the thumbnail task group? All good.
• Process thumbnail pixels sequentially instead of parallel? All good too.

🤯

The crash cause

Turns out, there’s nothing particularly wrong with the code above, it is “just” a surprising implementation detail of TBB task group cancellation, that is not intuitive at all.

It maybe makes sense if you would think really hard about “so, how would I actually implement task cancellation with nested parallelism?”, but most people do not think about this question every day.

What happens is this:

1. Each parallel_for creates a “task group context” (TBB task_group_context) as an on-stack / local variable. Our parallel_for_each(all_scene_geometries, build_geometry_bvh) above has just created one.
2. Whenever something “uses” a task group context, TBB “binds” it (whatever that is). As part of this “binding”, the context records the currently executing context as its parent (task_group_context.cpp:118), and adds itself into a per-thread list of live contexts (task_group_context.cpp:105). Turns out, Embree library also uses TBB internally, so building a BVH for a geometry does a bunch of parallel_for algorithms, and the TBB task_group_context objects get stored in our resulting BVH data. They all now point to the outer local context (created in previous step) as their “parent” though!
3. When using “Persistent Data” setting in Cycles, all our geometry BVH objects live for a long time; long after the initial “build all geometries in parallel” function has finished. They all have task_group_context objects that point to now-random stack data. So far, however, this does not cause any problems…
4. Later, we get to our thumbnails UI code. This has a completely unrelated task_group, that we call cancel() on. TBB then does cancellation propagation, and this walks a live-context list of every registered thread, and, for each context, walks its my_parent pointer chain up, looking for the context we are cancelling (task_group_context.cpp:200). This is where we find the still-alive task group created in step 2, and try to walk its parent into long-invalid memory since it was pointing to on-stack task group created in step 1. Note that this only happens if the task group that is cancelled itself has nested parallelism.

😤

I have a completely standalone repro demonstrating the issue here: https://github.com/aras-p/test_tbb_cancel

The lesson

Stop cancelling TBB tasks, I guess?

By now we know at least two things:

• In presence of task cancellation, parallel loops might not execute all of their iterations. And your code might need to be prepared to deal with that, even if you do not cancel any tasks! Someone way above in the call stack might.
• Long-lived task group contexts (as used by Embree) might record a surprise pointer to their caller parallel-for construct data, that can get stale. That pointer is only used when doing task cancellation. Even if cancellation happens on an entirely unrelated task group, in entirely unrelated code!

Luckily for Blender, it seems that task cancellation was used only in three places of the code. At least two of them used it for no good reason; I guess the authors saw “oh I can cancel a task group, sounds convenient” and used that. So another lesson would be, either do not expose task group cancellation functionality, or make the function sound much more scary than just cancel(), and have a giant comment warning the users of thing that lumbers slobberingly into sight and gropingly squeezes its geͣ͌͞l͓̻̥̓̿͐͐ă̜͆͜t̤ͥ͠i͕͚ͧͅṇ̽ͯ̀o̟u͚̯ͬ͐́͠s̷̜ͬ̓̏͆ g̴̼̹̼͊̅̓ŗ̶̝̟͌e̗͉̭ͤen͋ i̙̤̖͂͡_̷mme̷͈̺͙̓͞ͅn̥̾͐̃̒͊s̭iͩ̋ty̢͈ t̸̴̬͒͡ḥͬͧr̢̛ö̕u͒ğͯh̛͓̠̣̑ͣ͗ ṭ̂͂̑hͩë͕̻̙̐ b͋la̳̭̓c̤̓̓ͨk͖̠͉̽̃̚͠ d̘ͩoo̬̦͆ͣ̐͢r̻ͫͫ͛̚ẇ̧͕a̶̻͓͉y̡̬͍͊̏̓̂ i̠͑̈́n̘̪̋̓͛ͨt̨͓ͥ̋̀ͨ́o̡͆ͫ̿ t͓͆͛͌ͭͥhẽ͓ t͂ͧ̈ͮ́a̳̋̊ͣ̽ͭ͘ĩ͊n̶̿̓̏tͨͩe͓͋̓͡d͟ o̪̮̊̊u͑t̷̯̞̿̀͆͟šid̴͇͍̆̆͞ę̦͒͐͆ aͬ̍͆̊͋i̢͇̥̹̞̺ͥ̊͗ͫȓ̝̥̰̳͂̈́̚͢͡ ǫ̛̲̙̓́̀̇̈ͪ̕f̵̨͉ͥ͝ t͛̊ͨ̍h̪̪͋a̶̬͈͔̲̒̈́t p̡̨̗̗ͦͧ̾̎̚͠oͦ̍i̳̬͕̙̚s̥̜̗͑͜on̸̸̤̮ͨ̉͌̃ͫ c̢̡͇̳̲̻͐ͧȉ͕̔_t̯͙͉̑͟͞y ȏ̴̹̟̙ͭͧ͗f͕̠͆̎͐ m̪͈̞͙ͥ̿̽̅a̹͐ͫͭ̈́ͪ̒́d̸̶̞̣͔̤̟͚̼̔n̝̋e̫̩ͪ̎ͭ͟͡ss̴.

We will just remove usages of it: PR 160714, PR 160711 and PR 160709, and then un-expose TBB task group cancellation functionality from the rest of Blender codebase.

A tweet by @VehiclePhysics sparked my interest. It basically says:

For most math functions (Sqrt, Sin, Cos, Log, Pow…), prefer System.MathF over UnityEngine.Mathf. Unity’s Mathf casts to double, calls the double version, then converts back to float. System.MathF calls the float-native implementations directly. Less work, same result.

This advice is basically correct! But turns out, things are slightly more complicated.

Hidden double precision in Unity

The advice above applies to all UnityEngine.Mathf methods that deal with trigonometry (Sin, Cos, Tan, Asin, Acos, Atan, Atan2), exponentials (Sqrt, Pow, Exp, Log, Log10), rounding (Ceil, Floor, Round, CeilToInt, FloorToInt, RoundToInt), comparisons (Min, Max, Clamp, Clamp01) and others (Sign, SmoothStep, Gamma, Approximately, InverseLerp). About the only function it does not apply to is Mathf.Abs.

But… why? Well, because C#/.NET originally did not have single-precision methods for these sorts of math functions. The single precision System.MathF was introduced in .NET Core 2.0 (year 2017).

Now, you might have expected that almost ten years later, maybe Unity would have noticed this, and made them single precision? Alas, no. There could be potential backwards compatibility issues preventing that (or maybe not! see below).

You also might have guessed that Unity.Mathematics package, which was introduced (year 2019) as part of the whole DOTS push, and is modeled to be very similar to HLSL, would actually do single precision floating point for functions that look like single precision floating point… and that would be wrong too; for all the trigonometric and exponential functions like math.sqrt(float x) it routes that into the double precision C# implementation. Why? I don’t know.

But wait! There is way more double precision. The Mono C# runtime used in Unity does all math in double precision, everywhere. Yes, this means there is a ton of float⭤double conversions from in-memory representation to in-register representation, all over the place. I have first noticed this back in 2018 when doing a toy path tracer, and then Miguel de Icaza did an explanatory blog post, with plans outlined how to switch Mono to use actual floats for floats (yeah!).

“In Mono, decades ago, we made the mistake of performing all 32-bit float computations as 64-bit floats while still storing the data in 32-bit locations.”

Official Mono releases have switched to do that since then, but (I think) for backwards compatibility reasons Unity never enabled that functionality and kept everything at double precision so far.

Note however that the above only applies to Mono. The other two C# language/runtime implementations used across Unity today, IL2CPP and Burst, do not have the “everything is actually double precision” behavior. It is weird that Unity would not switch their Mono version to match; after all some of their main deployment platforms never use Mono (iOS, consoles, web)!

Let’s look at a square root

The above is fairly abstract, so let’s look at what actually happens with a very simple loop that sums up a bunch of square roots:

const int N = 10000000;
public static float UnityMathf(float v)
{
    for (int i = 0; i < N; ++i)
    {
        v += UnityEngine.Mathf.Sqrt(v); // classic Unity
        //v += System.MathF.Sqrt(v); // as advised by the tweet above
    }
    return v;
}

In the Unity editor (6000.0.76, but rough timings are the same on 2022.3, 6000.3 and 6000.6 versions), on Windows / Ryzen 5950X machine: UnityEngine.Mathf 282ms, System.MathF 186ms. Whoa indeed, this is way faster!

But hey! Back in 2018 we already found that Unity’s C# performance also very much depends on whether script debugging is enabled or not. Back then it was called “Editor Attaching” under preferences; these days it is this bad-contrast-in-light-theme Debug vs Release widget at lower right editor corner. In Release mode, in-editor timings are: UnityEngine.Mathf 242ms, System.MathF 149ms.

More square roots in more C# variants

To get a more complete picture, let’s also add a variant that uses the “new way of doing math” in Unity, i.e. the Unity.Mathematics package. And have timings for a player build that uses Mono, plus timings for an IL2CPP scripting backend. And while at it, also test performance of the same code under Burst compiler.

And for a complete picture, the same loop, using System.MathF.Sqrt (C#) or sqrtf() (C++) in non-Unity implementations / runtimes:

Summary of the above:

• 35 milliseconds to do this loop is “as good as it can get” on this machine, and that is achieved by C++ & .NET, and within Unity by using Burst + Unity.Mathematics, or when using IL2CPP, with either of Mathf.Sqrt or System.MathF.Sqrt. Under IL2CPP, there does seem to be some special code path that goes “oh this should actually be single precision square root” and generates underlying C++ code accordingly.
• System.MathF functions are not supported by Burst for some reason; if you try to use them you will get Burst compile errors. If you do not need Burst, then System.MathF is often faster. It does make it harder to move code to Burst though.
• Unity.Mathematics is often slightly better than the classic Mathf, except under IL2CPP, at least for the square root. IL2CPP does not seem to have special recognition of “oh this should be single precision square root” for it, and has other overheads too, see below.
• In the opposite behavior to IL2CPP, Burst does not seem to do “oh this should be single precision” for Mathf.Sqrt, but it does for Mathematics.math.sqrt at single precision.

Also fun fact? All the Unity implementations above print the result of the above loop as 24212990000000.0, which is curiously not a number that exists as a single precision float (closest floats that exist are 24212989280256.0 and 24212991377408.0). That’s one of the signs of “yeah some stuff is always doubles underneath, somewhere”. The non-Unity (C# .NET, C++) implementations print the result 24212987183104.0.

Welcome to the world! Things are never simple!

Code generation of the square root loops in detail

Mono, UnityEngine.Mathf.Sqrt

As the original tweet says, Unity’s Mathf.Sqrt is implemented like this: public static float Sqrt(float f) => (float)Math.Sqrt((double)f); – it just calls into double precision System.Math.Sqrt. But if you look at the actual JIT’ed machine code generated by Mono, you can see that there is way more float⭤double conversions going on.

I have used Sebastian Schöner’s Asm Explorer tool to see the generated code. Given this C# code:

const int N = 10000000;
public static float UnityMathf(float v)
{
    for (int i = 0; i < N; ++i)
    {
        v += UnityEngine.Mathf.Sqrt(v);
    }
    return v;
}

the loop body ends up being this:

loop:
movss xmm0, dword [rsp+0x10]     ; xmm0 = v, as float
cvtss2sd xmm0, xmm0              ; xmm0 = (double)v, left side of v + sqrt(v)

movss xmm1, dword [rsp+0x10]     ; xmm1 = v, as float again, argument for sqrt
cvtss2sd xmm1, xmm1              ; xmm1 = (double)v

cvtsd2ss xmm5, xmm1              ; xmm5 = (float)(double)v, rounded back to float
movss [rsp+0x8], xmm5            ; store temporary float argument

movss xmm1, dword [rsp+0x8]      ; xmm1 = temporary float argument
cvtss2sd xmm1, xmm1              ; xmm1 = (double)temporary float

movsd [rsp-0x8], xmm1            ; store double for x87 sqrt input
fld qword [rsp-0x8]              ; push double onto x87 stack
fsqrt                            ; ST(0) = sqrt(ST(0))
fstp qword [rsp-0x8]             ; store sqrt result as double and pop x87 stack

movsd xmm1, qword [rsp-0x8]      ; xmm1 = sqrt result, as double
cvtsd2ss xmm1, xmm1              ; xmm1 = (float)sqrt result
cvtss2sd xmm1, xmm1              ; xmm1 = (double)(float)sqrt result

cvtsd2ss xmm5, xmm1              ; xmm5 = sqrt result rounded to float
movss [rsp+0x8], xmm5            ; store temporary sqrt float

movss xmm1, dword [rsp+0x8]      ; xmm1 = temporary sqrt float
cvtss2sd xmm1, xmm1              ; xmm1 = (double)temporary sqrt float

cvtsd2ss xmm5, xmm1              ; xmm5 = sqrt result rounded to float again
movss [rsp+0x8], xmm5            ; store temporary sqrt float again

movss xmm1, dword [rsp+0x8]      ; xmm1 = temporary sqrt float again
cvtss2sd xmm1, xmm1              ; xmm1 = (double)sqrt result

addsd xmm0, xmm1                 ; xmm0 = (double)v + (double)Mathf.Sqrt(v)

cvtsd2ss xmm5, xmm0              ; xmm5 = final iteration result, round to float
movss [rsp+0x10], xmm5           ; v = final iteration result

inc esi                          ; ++i
cmp esi, 0x989680                ; compare i against 10000000
jl loop                          ; if i < N, continue loop

If this were C#, it would be like v += UnityEngine.Mathf.Sqrt(v) actually expands to:

double lhs = (double)v;

double t0 = (double)v;
float t1 = (float)t0;
float stackFloat0 = t1;
float t2 = stackFloat0;
double sqrtInput = (double)t2;

double stackDouble0 = sqrtInput;
double sqrtDouble = X87_Fsqrt(stackDouble0); // represents x87 FPU fsqrt instruction

float t3 = (float)sqrtDouble;
double t4 = (double)t3;
float t5 = (float)t4;
float stackFloat1 = t5;
float t6 = stackFloat1;
double t7 = (double)t6;
float t8 = (float)t7;
float stackFloat2 = t8;
float t9 = stackFloat2;
double rhs = (double)t9;

double sum = lhs + rhs;
float result = (float)sum;
v = result;

That’s… not exactly great, to put it mildly. Unity is planning to switch to “actual .NET” (CoreCLR) really soon now (see Path to CoreCLR GDC 2026 talk) and codegen should get much better then. Meanwhile, I am rediscovering the same things as what Sebastian Schöner did, but he is also trying to do something about it – see Better codegen for Unity games on Mono blog post.

Using Unity.Mathematics.math.sqrt is a tiny bit better codegen than above, but not by much.

Mono, System.MathF.Sqrt

const int N = 10000000;
public static float UnityMathf(float v)
{
    for (int i = 0; i < N; ++i)
    {
        v += System.MathF.Sqrt(v);
    }
    return v;
}

the loop body ends up being this:

loop:
movss xmm0, dword [rbp-0x10]     ; xmm0 = v, as float
cvtss2sd xmm0, xmm0              ; xmm0 = (double)v
movsd [rbp-0x18], xmm0           ; save old v as double for later addition

movss xmm0, dword [rbp-0x10]     ; xmm0 = v, as float again, argument for MathF.Sqrt
cvtss2sd xmm0, xmm0              ; xmm0 = (double)v

cvtsd2ss xmm0, xmm0              ; xmm0 = (float)(double)v, argument to MathF.Sqrt
nop                              ; padding / alignment / patchpoint artifact

mov r11, 0x22494ee3918           ; r11 = JIT trampoline address for System.MathF.Sqrt(float)
call r11                         ; call MathF.Sqrt(float), argument in xmm0, return float in xmm0

cvtss2sd xmm1, xmm0              ; xmm1 = (double)MathF.Sqrt(v)

movsd xmm0, qword [rbp-0x18]     ; xmm0 = saved old v, as double
addsd xmm0, xmm1                 ; xmm0 = (double)old_v + (double)sqrt_v

cvtsd2ss xmm5, xmm0              ; xmm5 = final iteration result rounded to float
movss [rbp-0x10], xmm5           ; v = final iteration result

inc esi                          ; ++i
cmp esi, 0x989680                ; compare i against 10,000,000
jl loop                          ; if i < N, continue loop

and the assembly of the actual System.MathF.Sqrt function is:

xorps xmm1, xmm1                 ; xmm1 = 0.0f
ucomiss xmm1, xmm0               ; compare 0.0f with input
ja handlefail                    ; if 0.0f > input, input is negative: go handle failure/NaN path
sqrtss xmm0, xmm0                ; xmm0 = sqrtss(xmm0), scalar single-precision sqrt
ret                              ; return sqrt result in xmm0

handlefail:
; some code that handles failures/NaNs

it is effectively this:

static float MathF_Sqrt_Call(float x)
{
    if (0.0f > x)
        return MathF_Sqrt_SlowPath(x);
    return Sse_SqrtScalarSingle(x); // sqrtss instruction
}

// ...
double lhs = (double)v;

double t0 = (double)v;
float sqrtArg = (float)t0;
float sqrtResult = MathF_Sqrt_Call(sqrtArg);
double rhs = (double)sqrtResult;

double sum = lhs + rhs;
float result = (float)sum;
v = result;

There are still a bunch of float⭤double conversions! But way fewer, and instead of using the ancient x87 FPU, this now uses the scalar SSE square root instruction.

Burst, UnityEngine.Mathf.Sqrt

Under Burst, the v += UnityEngine.Mathf.Sqrt(v) inner loop part faithfully translates to:

vcvtss2sd   xmm1, xmm0, xmm0 ; convert float→double
vsqrtsd     xmm1, xmm1, xmm1 ; scalar double precision square root
vcvtsd2ss   xmm1, xmm1, xmm1 ; convert double→float
vaddss      xmm0, xmm0, xmm1 ; float +=

i.e. it does pretty much what you would expect, given Mathf.Sqrt implementation.

Burst, Unity.Mathematics.math.sqrt

The v += Unity.Mathematics.math.sqrt(v) under Burst translates to just:

vsqrtss     xmm1, xmm0, xmm0 ; scalar single precision square root
vaddss      xmm0, xmm0, xmm1 ; float +=

This is basically what you would want to happen.

This is somewhat curious though, since underlying math.sqrt code is actually public static float sqrt(float x) { return (float)System.Math.Sqrt((float)x); } – i.e. without Burst, it does end up calling into double precision function. But Burst gives this some sort of special treatment, that it does not do for the previous case, I guess.

And again, no System.MathF.Sqrt test with Burst, since it just fails if you try to use that.

IL2CPP, UnityEngine.Mathf.Sqrt

Unity’s IL2CPP scripting backend translates .NET bytecode into C++, and then relies on a regular C++ compiler to carry out optimizations.

For the Mathf.Sqrt code path, it does seem to actually give it special treatment – it does not call the double precision square root, even if on C# level it does do double precision. This is the opposite of what Burst does, and I guess this is another example of “you ship your org chart” in action.

The inner loop in generated C++ code is:

float L_0 = ___0_v;
float L_1 = ___0_v;
float L_2;
L_2 = sqrtf(L_1);
___0_v = (float)il2cpp_codegen_add(L_0, L_2); // template function, just + for simple types

which then the C++ compiler (MSVC 2022 v17.14, Release build config) actually unrolls to do ten square roots per iteration, with each square root snippet being this:

xorps       xmm1, xmm1     ; xmm1 = 0.0f
ucomiss     xmm1, xmm6     ; compare 0.0f with v
ja          edgecase       ; if 0.0f > v: use sqrtf fallback function
xorps       xmm0, xmm0     ; xmm0 = 0.0f
sqrtss      xmm0, xmm6     ; xmm0 = sqrt(v), scalar single-precision sqrt
jmp         end
edgecase:
movaps      xmm0, xmm6     ; xmm0 = v, argument for sqrtf
call        sqrtf          ; call C runtime sqrtf
end:
addss       xmm6, xmm0     ; v += sqrtResult

This is not a simple “just use sqrtss”, it only uses the instruction for valid inputs, and calls into “full” function for others (to set errno or deal with exceptions, I guess). You could argue that this is less optimal codegen than what Burst does, in practice on this benchmark it does not matter though.

IL2CPP, System.MathF.Sqrt

Now, for System.MathF the IL2CPP codegen is slightly different:

il2cpp_codegen_runtime_class_init_inline(MathF_longGUID_il2cpp_TypeInfo_var);
float L_0 = ___0_v;
float L_1 = ___0_v;
float L_2;
L_2 = sqrtf(L_1);
___0_v = (float)il2cpp_codegen_add(L_0, L_2); // template function, just + for simple types

– why yes, that is the il2cpp_codegen_runtime_class_init_inline call inside the hot inner loop. What that does, is it checks some flag and if it is not set, calls some other function. Some sort of “lazy C# class initialization”, that for some reason is not needed in the previous case, but is needed here.

In assembly, this looks very much like above, except now the loop body is not “tiny enough” so MSVC compiler does not do ten square roots per each actual loop iteration; it does only one. And before each square root, it does this:

mov         rcx,qword ptr [MathF_longGUID_il2cpp_TypeInfo_var]  
cmp         dword ptr [rcx+0E4h],0  
jne         inited
call        il2cpp_codegen_runtime_class_init
inited:

Now again, for this particular benchmark it does not matter (the memory address it checks is very much in the cache, and the branch is perfectly predictable). But if you are calling System.MathF.Sqrt outside of tiny inner loops, then each.and.every.call will have this extra memory fetch and a branch.

IL2CPP, Unity.Mathematics.math.sqrt

For the Mathematics.math.sqrt case, things get slightly weirder under IL2CPP: 1) instead of one “some sort of lazy initialization” branch like in case above, now it has two branches for each and every call, and 2) the actual square root is done in double precision.

Generated C++ code:

IL2CPP_MANAGED_FORCE_INLINE IL2CPP_METHOD_ATTR float math_sqrt_longGUID_inline (float x, const RuntimeMethod* method) 
{
  static bool s_Il2CppMethodInitialized;
  if (!s_Il2CppMethodInitialized)
  {
    il2cpp_codegen_initialize_runtime_metadata((uintptr_t*)&Math_longGUID_il2cpp_TypeInfo_var);
    s_Il2CppMethodInitialized = true;
  }
  {
    il2cpp_codegen_runtime_class_init_inline(Math_longGUID_il2cpp_TypeInfo_var);
    double l1 = sqrt((double)x);
    return (float)l1;
  }
}

which then translates into this assembly for the inner loop:

cmp         byte ptr [s_Il2CppMethodInitialized],0  
jne         inited1
lea         rcx,[Math_longGUID_il2cpp_TypeInfo_var]  
call        il2cpp_codegen_initialize_runtime_metadata
mov         byte ptr [s_Il2CppMethodInitialized],1  
inited1:
mov         rcx,qword ptr [Math_longGUID_il2cpp_TypeInfo_var]  
cmp         dword ptr [rcx+0E4h],0  
jne         inited2
call        il2cpp_codegen_runtime_class_init
inited2:
xorps       xmm1,xmm1  
xorps       xmm0,xmm0  
cvtss2sd    xmm1,xmm6  
ucomisd     xmm0,xmm1  
ja          edge_case
sqrtpd      xmm0,xmm1  
jmp         iter_end
edge_case:
movaps      xmm0,xmm1  
call        sqrt
iter_end:
cvtsd2ss    xmm0,xmm0  
addss       xmm6,xmm0 

Again, for this benchmark the two extra branches do not matter, but they might if you are calling math.sqrt not from inside of a tiny loop body. What does matter, and why under IL2CPP this is slower, is that the square root is done at double precision.

So there! Unity math is complex!

Well, that was something. Is the original advice of prefer System.MathF over UnityEngine.Mathf valid? Yes, unless you want Burst; there it simply does not work.

My takeaways:

• I hope the upcoming switch to .NET / CoreCLR will clear up a lot of that mess, especially in the “even if you don’t spell out doubles anywhere in your code, Mono does everything in doubles in Unity”. And even without double precision, the Mono codegen is… not great.
• Unity is quite inconsistent in how it treats precision of various math functions. Some of them are implemented as-if they were double precision, but IL2CPP and Burst magically treat them as single precision. Sometimes IL2CPP and Burst disagree on which ones get the special treatment.
  • Given that CoreCLR switch will have some potential backwards compat breakages anyway, I hope Unity will sanitize the math functions precision treatments in the same go.
• It would be nice if you could use “functions that look & feel the same” (like UnityEngine.Mathf.Sqrt, System.MathF.Sqrt and Unity.Mathematics.math.sqrt) as being exactly equivalent, with no preferential treatment of one vs. the other. That is very much not the case today however, and what’s worse, there is no single answer for “which one is best”. It all depends whether you use IL2CPP or Burst, or both, or neither!
• If you want best performance now, use Burst and Mathematics maths.
• Also, you might want to look into Sebastian’s cpp2better, that is aimed at improving IL2CPP codegen. I have not evaluated it in this post however.