qznc 
Posted in reply to bearophile
| On Tuesday, 3 May 2011 at 20:51:37 UTC, bearophile wrote:
> Sean Cavanaugh:
>
>> In many ways the biggest thing I use regularly in game development that I would lose by moving to D would be good built-in SIMD support.
>
> Don has given a nice answer about how D2 plans to face this.
>
> To focus more what Don was saying I think a small exaple will help. This is a C implementation of one Computer Shootout benchmarks, that generates a binary PPM image of the Mandelbrot set:
>
> http://shootout.alioth.debian.org/u32/program.php?test=mandelbrot&lang=gcc&id=4
>
> This is an important part of that C version:
>
>
> typedef double v2df __attribute__ ((vector_size(16))); /* vector of two doubles */
> const v2df zero = { 0.0, 0.0 };
> const v2df four = { 4.0, 4.0 };
>
> // Constant throughout the program, value depends on N
> int bytes_per_row;
> double inverse_w;
> double inverse_h;
>
> // Program argument: height and width of the image
> int N;
>
> // Lookup table for initial real-axis value
> v2df *Crvs;
>
> // Mandelbrot bitmap
> uint8_t *bitmap;
>
> static void calc_row(int y) {
> uint8_t *row_bitmap = bitmap + (bytes_per_row * y);
> int x;
> const v2df Civ_init = { y*inverse_h-1.0, y*inverse_h-1.0 };
>
> for (x = 0; x < N; x += 2) {
> v2df Crv = Crvs[x >> 1];
> v2df Civ = Civ_init;
> v2df Zrv = zero;
> v2df Ziv = zero;
> v2df Trv = zero;
> v2df Tiv = zero;
> int i = 50;
> int two_pixels;
> v2df is_still_bounded;
>
> do {
> Ziv = (Zrv * Ziv) + (Zrv * Ziv) + Civ;
> Zrv = Trv - Tiv + Crv;
> Trv = Zrv * Zrv;
> Tiv = Ziv * Ziv;
>
> // All bits will be set to 1 if 'Trv + Tiv' is less than 4
> // and all bits will be set to 0 otherwise. Two elements
> // are calculated in parallel here.
> is_still_bounded = __builtin_ia32_cmplepd(Trv + Tiv, four);
>
> // Move the sign-bit of the low element to bit 0, move the
> // sign-bit of the high element to bit 1. The result is
> // that the pixel will be set if the calculation was
> // bounded.
> two_pixels = __builtin_ia32_movmskpd(is_still_bounded);
> } while (--i > 0 && two_pixels);
>
> // The pixel bits must be in the most and second most
> // significant position
> two_pixels <<= 6;
>
> // Add the two pixels to the bitmap, all bits are
> // initially zero since the area was allocated with calloc()
> row_bitmap[x >> 3] |= (uint8_t) (two_pixels >> (x & 7));
> }
> }
>
>
> GCC 4.6 compiles the inner do-while loop of calc_row() to just this very clean assembly, that in my opinion is quite _beautiful_, it shows one of the most important final purposes of a good compiler:
>
> L9:
> subl $1, %ecx
> addpd %xmm0, %xmm0
> mulpd %xmm0, %xmm1
> movapd %xmm4, %xmm0
> addpd %xmm6, %xmm1
> addpd %xmm5, %xmm0
> subpd %xmm3, %xmm0
> movapd %xmm1, %xmm3
> movapd %xmm0, %xmm4
> mulpd %xmm1, %xmm3
> mulpd %xmm0, %xmm4
> movapd %xmm3, %xmm2
> addpd %xmm4, %xmm2
> cmplepd %xmm7, %xmm2
> movmskpd %xmm2, %ebx
> je L18
> testl %ebx, %ebx
> jne L9
>
>
> Those addpd, subpd, mulpd, movapd, etc, instructions work on pairs of doubles (those v2df). And the code uses the cmplepd and movmskpd instructions too, in a very clean way, that I think not even GCC 4.6 is normally able to use by itself. A good language + compiler have many purposes, but producing ASM code like that is one of the most important purposes, expecially if you write numerical code.
>
> A numerical programmer really wants to write code that somehow produces equally clean and powerful code (or better, using AVX 256-bit registers and 3-way instructions) in numerical processing kernels (often such kernels are small, often just bodies of inner loops).
>
> D2 allows to write code almost as clean as this C one (but I think currently no D compiler is able to turn this into clean inlined addpd, subpd, mulpd, movapd instructions. This is a compiler issue, not a language one):
>
> v2df Zrv = zero;
> ...
> Ziv = (Zrv * Ziv) + (Zrv * Ziv) + Civ;
> Zrv = Trv - Tiv + Crv;
> Trv = Zrv * Zrv;
> Tiv = Ziv * Ziv;
>
>
> In D it becomes:
>
> double[2] Zrv = zero;
> ...
> Ziv[] = (Zrv[] * Ziv[]) + (Zrv[] * Ziv[]) + Civ[];
> Zrv[] = Trv[] - Tiv[] + Crv[];
> Trv[] = Zrv[] * Zrv[];
> Tiv[] = Ziv[] * Ziv[];
>
>
> But then how do you write this in a clean way in D2/D3?
>
> do {
> ...
> is_still_bounded = __builtin_ia32_cmplepd(Trv + Tiv, four);
> two_pixels = __builtin_ia32_movmskpd(is_still_bounded);
> } while (--i > 0 && two_pixels);
>
>
>
> Using those __builtin_ia32_cmplepd() and __builtin_ia32_movmskpd() is not easy, so there is a tradeoff between allowing easy to write code, and giving power. So it's acceptable for a language to give a bit less power if the code is simpler to write. Yet, in a system language if you don't give people a way to produce ASM code as clean as the one I've shown in the inner loops of numerical processing code, some D2 programmers will be forced to write down inline asm, and that's sometimes worse than using intrinsics like __builtin_ia32_cmplepd().
>
> Writing efficient inner loops is very important for numerical processing code, and I think numerical processing code is important for D2.
>
> Time ago I have suggested to extend the D2 vector operations to code like this, but I think this is not enough still:
>
> float[4] a, b, c, d;
> c = a[] == b[];
> d = a[] >= b[];
>
> Bye,
> bearophile
Just found this old post, since I'm tuning mandelbrot.d right now [0].
The good news: LDC produces code, which is quite close to the C version.
mulsd %xmm6,%xmm4
subsd %xmm1,%xmm7
addsd %xmm4,%xmm4
addsd %xmm5,%xmm7
addsd %xmm0,%xmm4
movaps %xmm7,%xmm6
mulsd %xmm6,%xmm6
movaps %xmm4,%xmm2
mulsd %xmm2,%xmm2
movaps %xmm2,%xmm1
addsd %xmm6,%xmm1
ucomisd %xmm1,%xmm3
jb 4026f0 <_D10mandelbrot11computeLineFNaNbNfmiZAa+0x130>
jl 402680 <_D10mandelbrot11computeLineFNaNbNfmiZAa+0xc0>
Even better, the code is produce from the following (inlined!) source,
which is pretty much the mathematical definition.
for(auto i = 0; i < iter && norm(Z) <= lim; i++)
Z = Z*Z + C;
The bad news: cmplepd and movmskpd are not used. Is that possible somehow four years later?
The gcc code is roughly twice as fast at the moment, but I don't know if cmplepd and movmskpd is the only thing missing.
[0] https://github.com/qznc/d-shootout
|
Permalink
Reply
