Optimizing the Particle Life: From 400 to 4M particles

Post by Twirrim » Mon, Mar 18, 2024, 8:26 PM UTC

bscottmay wrote: post #39

Twirrim wrote: post #30

The thing that strikes me as particularly weird about this is the use of numpy there. In other languages, they're using native code. In python they're reaching out to numpy, which is a great library, but not awesome inside a hot loop unless you're keeping the operation you're carrying out within numpy itself. This means, right in that hot loop, they're doing a lot of translating of numbers between python representati…

Yeah, an all-numpy version runs in less than 1ms on my M1 air import numpy as np l = 10_000 t = np.empty(l, dtype=np.float32) j = np.arange(l) t = 0.02 \* j t *= (0.03 * j) t -= (0.04 \* j) t /= 0.05 \* (j + 1)

Post by jasonjmcghee » Mon, Mar 18, 2024, 8:29 PM UTC

wzdd wrote: post #6

hyperbrainer wrote: post #2

Why was C slower though?

Simply compiling it with -O3 produces something which completes in half the time of the JavaScript version (350ms for C, 750ms for JS), so perhaps that. Edit for Twirrim: on this system (Ryzen 7, gcc 11): "-O3": 350ms; "-O3 -march=native": 208ms; "-O2": 998ms; "-O2 -march=native": 1040ms. Edit 2: Interestingly, changing the C from float to double produces a 3.5x speedup, taking the time elapsed (with "-O3 -march=nati…

Post by [deleted] » Mon, Mar 18, 2024, 8:39 PM UTC

a1369209993 wrote: post #45

abainbridge wrote: post #40

Earlier quoted context omitted.

One optimization for the C code is to put "f" suffixes on the floating point constants. For example convert this line: t[i] += 0.02 * (float)j; to: t[i] += 0.02f * (float)j; I believe this helps because 0.02 is a double and doing double * float and then converting the result to float can produce a different answer to just doing float * float. The compiler has to do the slow version because that's what you asked for.…

> I believe this helps because 0.02 is a double and [...] can produce a different answer In principle, not quite. The real/unavoidable(-by-the-compiler) problem is that 0.02 is a not a diadic rational (not representable exactly as some integer over a power of two). So its representation (rounded to 52 bits) as a double is a different real number than its representation (rounded to 23 bits) as a float. (This is the sa…

Post by abainbridge » Mon, Mar 18, 2024, 8:43 PM UTC

a1369209993 wrote: post #45

abainbridge wrote: post #40

Earlier quoted context omitted.

One optimization for the C code is to put "f" suffixes on the floating point constants. For example convert this line: t[i] += 0.02 * (float)j; to: t[i] += 0.02f * (float)j; I believe this helps because 0.02 is a double and doing double * float and then converting the result to float can produce a different answer to just doing float * float. The compiler has to do the slow version because that's what you asked for.…

> I believe this helps because 0.02 is a double and [...] can produce a different answer In principle, not quite. The real/unavoidable(-by-the-compiler) problem is that 0.02 is a not a diadic rational (not representable exactly as some integer over a power of two). So its representation (rounded to 52 bits) as a double is a different real number than its representation (rounded to 23 bits) as a float. (This is the sa…

    #include 
    #include 

    float mul_as_float(float t) {
      t += 0.02f * (float)17;
      return t;
    }

    float mul_as_double(float t) {
      t += (double)0.02f * (float)17;
      return t;
    }

    int main() {
        while (1) {
            unsigned r = rand();
            float t = *((float*)&r);

            float result1 = mul_as_float(t);
            float result2 = mul_as_double(t);
            if (result1 != result2) {
                printf("Counter example when t is %f (0x%x)\n", t, *((unsigned*)&t));
                printf("result1 is %f (0x%x)\n", result1, *((unsigned*)&result1));
                printf("result2 is %f (0x%x)\n", result2, *((unsigned*)&result2));
                return 0;
            }
        }
    }

    Counter example when t is 0.000000 (0x3477d43f)
    result1 is 0.340000 (0x3eae1483)
    result2 is 0.340000 (0x3eae1482)

Post by abainbridge » Mon, Mar 18, 2024, 9:35 PM UTC

a1369209993 wrote: post #58

abainbridge wrote: post #57

Earlier quoted context omitted.

OK, right, that's the clarity of thought I was missing. But in this case the compiler still misses the optimization with '(double)0.02f'. https://godbolt.org/z/az7819nKM I think this is because the optimization isn't safe. I wrote a program to find a counter example to your claim that "the optimization should in theory still apply". It found one. Here's the code: #include #include float mul_as_float(float t) { t += 0…

On my machine, the complier constant-folds the multiplication, producing a single-precision add for `mul_as_float` and a convert-t-to-double, double-precision-add, convert-sum-to-single for `mul_as_double`. I missed the `+=` in your original comment, but adding a float to a double does implicitly promote it like that, so you'd actually need: t += (float)((double)0.02f * (float)17); to achieve the "and then converting…