Please note I used the following code to test the exp and log functions:

void test_ext(uniform float x, uniform bool verbose)
{
    const uniform float res = log(x);
    const uniform float dres = (uniform float)log((uniform double)x);
    const uniform float inf = floatbits(0x7F800000);

    if(verbose)
        print("%: %\n", x, res);

    if(isnan(res) != isnan(dres))
        print("ERROR: NaN mismatch: x=%\n", x);
    else if(res != dres && (abs(res) == inf || abs(dres) == inf))
        print("ERROR: Inf mismatch: x=% res=% dres=%\n", x, res, dres);
    else if(res != dres && abs(res-dres)/dres >= 130e-7) // Usual error check: >= 4~5 ULP
        print("ERROR: value mismatch: x=% res=% dres=% error_x1000=%\n", x, res, dres, abs(res-dres)*1000/dres);

    if(verbose)
        print("\n");
}

void test(uniform float x)
{
    test_ext(x, true);
}

extern "C" uniform int main()
{
    uniform bool full_check = true;
    uniform bool measure_perf = true;

    uniform float pinf  = floatbits(0x7F800000);
    uniform float ninf  = floatbits(0xFF800000);
    uniform float pnan  = floatbits(0x7FC00000);
    uniform float nnan  = floatbits(0xFFC00000);
    uniform float pnan2 = floatbits(0x7FC00123);

    test(0.0f);
    test(-0.0f);
    test(1.0f);
    test(1.0000005f); // very imprecise for initial log
    test(1.0000010f);
    test(1.0000011f); // very imprecise for initial log
    test(1.0000012f);
    test(-1.0f);
    test(-10.0f);
    test(10.0f);
    test(1.15e+5);
    test(1.15e+36);
    test(1.15e+42);
    test(1.15e-36);
    test(pinf);
    test(ninf);
    test(pnan);
    test(nnan);
    test(pnan2);

    if(full_check)
    {
        uniform int64 maxi = ((uniform int64)1) << 32;
        for(uniform int64 i=0; i<maxi; ++i)
        {
            test_ext(floatbits((uniform int32)i), false);

            if(i % 100000000 == 0)
                print("%/100\n", i*100/maxi);
        }
    }

    if(measure_perf)
    {
        uniform int32 maxi = (((uniform int64)1) << 31) - 1; // skip negative numbers and few NaN to avoid overflows
        float checksum = 0.0;
        const uniform int64 startTime = clock();

        foreach(i=0...maxi)
            checksum += log(floatbits(i));

        const uniform int64 endTime = clock();
        print("checksum: %, clock: %e9\n", reduce_add(checksum), (uniform float)(endTime-startTime)*1e-9);
    }
}

It may be useful for another use-case (e.g. other math functions) or possibly for tests. It is quite convenient to check all possible float values (since tests based on random numbers are not really reliable).

Thank you for you contribution!

I suggest splitting this PR into two because it addresses two different problems.

To make it look neat, I suggest formatting the code (see the output of this job) and formatting the git message as follows: the first line is a short description, e.g., stdlib: fix overflow in exp(float) then an empty line, followed by a detailed description of any length, formatted to a width of 80 characters.

Could you wrap your test code into a trivial microbenchmark, similar to what we have here: ISPC Benchmarks? This will allow us to retain and integrate your testing code effectively.

Does your change anyhow utilizes Sleef algorithms?

Could you wrap your test code into a trivial microbenchmark, similar to what we have here: ISPC Benchmarks? This will allow us to retain and integrate your testing code effectively.

We already have it for pow, log, and exp:

Also I measured the impact of the changes on Apple M3 chip with neon-i32x4 target. exp is not affected. log is slower approximately 50%.

Does your change anyhow utilizes Sleef algorithms?

No. This is my own modification.

Could you wrap your test code into a trivial microbenchmark, similar to what we have here: ISPC Benchmarks? This will allow us to retain and integrate your testing code effectively.

We already have it for pow, log, and exp

Ok, so there is finally nothing to add in the microbenchmark, isn't it?

Ok, so there is finally nothing to add in the microbenchmark, isn't it?

Yep, no changes.

I improved the speed of the log function. It is now "only" 25% slower than the initial code for the target avx2-i32x8 and 55% for the target avx2-i32x16. This mitigate the strange issue related to the Inf (the assembly code seems minor but the performance impact is surprisingly huge). I am going to investigate if the avx2-i32x16 target can be optimized further. I am not sure there is much more to do now regarding performance of the target avx2-i32x8 (or Neon). I think 25% is the price to pay for the fixes.

I tried to understand what was more precisely the numerical error so to use a more efficient approach for the correction close to x=1 and it worked well. The new approach is not only faster for the avx2-i32x8 target but also significantly more accurate. I described the approach more precisely in comments of the submitted ISPC file.

In the end, the new code is only 8% slower than the initial one for the avx2-i32x8 target and still 50% slower with the avx2-i32x16 (due to register spilling and apparently LLVM optimization issues too). The new code is slightly less precise (only few ULP) for most values, but much more precise for values close to x=1 : about 10-15 times better! Overall the new code has a relative precision of <10 ULP in practice (so 1_000-10_000 time more than the initial code).

I also fixed the formatting issues.
I am going to certainly remove the exp fix from this PR and open a new one as proposed by @nurmukhametov .

It should be fine now and ready to be merged :) !

In the end, the new code is only 8% slower than the initial one for the avx2-i32x8 target and still 50% slower with the avx2-i32x16 (due to register spilling and apparently LLVM optimization issues too).

It looks like if we move exceptional_result definition up like this:

@@ -3477,6 +3477,7 @@ __declspec(safe) static inline float log(float x_full) {
         const bool use_nan = !(x_full >= 0.);
         const bool use_inf = x_full == 0. || x_full == inf;
         const bool exceptional = use_nan || use_inf;
+        const float exceptional_result = select(use_nan, NaN, select(x_full == inf, x_full, -inf));
         const float one = 1.0;
 
         const float patched = select(exceptional, one, x_full);
@@ -3509,7 +3510,6 @@ __declspec(safe) static inline float log(float x_full) {
 
         result = log_from_exponent - x * result;
 
-        const float exceptional_result = select(use_nan, NaN, select(x_full == inf, x_full, -inf));
         return select(exceptional, exceptional_result, result);
     }
 }

then no spills are generated. Could you try this?

In practice, there are 6 spills in the current code and 7 spill if I move the line. The x16 target code is also 1% slower (not really significant). The x8 version seems not affected (as expected). Thus, it unfortunately does not help on my machine.

By the way, note that I saw different results between Godbolt and my local ISPC (with no other change than this PR) and I do not know why yet. Did you use that to check spilling?

I expected LLVM to track the dependency and perform a good register allocation so to reduce spilling as much as possible. In practice, such problem are not rare in ISPC with double-pumping. That being said, here, it is difficult for LLVM to produce a good code (except on AVX-512 thanks to cheap broadcasting embedded in instructions) based on the high number of registers to load and the number of variables. Still, there is a room for improvement since the x16 target is significantly slower than the x8 target and I think LLVM could theoretically just reorder most instructions so to generate operation on the first x8 items and then the second group. This is mostly the case but not totally. However, in practice, the boolean operation results in fused instructions preventing partially such optimization. This is a bit complicated since I think merging can be better when people do a lot of boolean operations on SSE/AVX-1/AVX-2/Neon targets while it is not a good idea otherwise except maybe in pathological cases (maybe it can reduce register pressure). I am not sure ISPC can make good decision while LLVM can certainly (thanks to a global optimization step). I am clearly not an expert of this so I may completely miss-understand how such thing work or how things should be improved. I think this problem is related to this issue. Still, I am not sure this is the only issue here.

For example, generated instructions are sub-optimal when spilling happens:

vblendvps	ymm2, ymm2, ymm6, ymm5
vblendvps	ymm2, ymm14, ymm2, ymm4
vmovups	ymm3, ymmword ptr [rsp - 112]   # 32-byte Reload
vaddps	ymm3, ymm3, ymm1
vmovups	ymmword ptr [rsp - 112], ymm3   # 32-byte Spill
vmovups	ymm1, ymmword ptr [rsp - 112]   # 32-byte Reload    <--------------------
vaddps	ymm0, ymm0, ymm2
add	eax, 16
cmp	eax, 2147483632
jne	.LBB109_1

Here we can see a useless reload which can be replaced by a basic vmovps ymm1, ymm3.

I think spilling can can responsible for up to a 30% slowdown on my machine so it may not fully explain the performance drop of the double-pumping version...

In practice, there are 6 spills in the current code and 7 spill if I move the line. The x16 target code is also 1% slower (not really significant). The x8 version seems not affected (as expected). Thus, it unfortunately does not help on my machine.

By the way, note that I saw different results between Godbolt and my local ISPC (with no other change than this PR) and I do not know why yet. Did you use that to check spilling?

I used my local build but the generated code looks similar to compiler-explorer's code link. I suggest to check which LLVM version you use. Release ISPC and trunk ISPC in compiler-explorer are built with the patched LLVM (see llvm_patches dir in repo). But I don't see difference on my machine between patched LLVM 17 and not-patched LLVM 18 for this example.

TL;DR: the overhead of the target avx2-i32x16 (on Zen2) appears to be mainly the price to pay for the much better precision and correct support of special value. I do not think we can reduce it much more (certainly not more than 10%) unless a better approach is found (and I do not think there is a much better approach both reaching a <5 ULP precision and supporting correctly special values).

I profiled more intensively the double pumping implementation to better know what is happening on Zen2 with this PR. The code is bound by the number of instruction/uops which is about 50% higher with the PR. Since the IPC is good both before/after the modification (~2 FP uops/cycle), this makes sense to get a 50% slowdown.

Here are few difference between the two versions for the target avx2-i32x16:

• the number of vblendvps instruction increases from 6 to 12 : this is normal because of the patch close to 0 so to increase the precision. Indeed, 3 new select are done (2 for the patch close to zero and 1 for the support of Inf).
• 2 FMA instructions seems to be split in sub+mul instruction (I do not know why and if they can theoretically be fused).
• the number of comparison instructions increases from 4 to 8 (new vcmpnleps/vcmpnltps) : this is normal due to the better support of Inf (first additional check) and the check of values close to 0 (second additional check).
• the number of vxorps+vpor increases from 4 to 8 overall : this is certainly due to the check too though the vxorps could certainly be optimized out.
• the number of mov-like instructions increases from 14 to 19 : this is due to register spilling (3 spill vs 6 spill) which is due to a higher register pressure itself certainly due to more constant to load and the additional operations. There is maybe a way to reduce it further but certainly not to reach the same amount as before.
• 2 additional vandps instructions are generated : this is due to the check close to 0 (i.e. the call to abs).
• few additional new instructions (x7) like vpshufd, vpackssdw and vpacksswb and vextracti128 are generated and I think this is missed optimization of LLVM/ISPC related to boolean operations. Indeed, we pack boolean of YMM registers so to perform the logical-or and then unpack the XMM registers so to get YMM ones back while we could just perform YMM operations when the number of booleans operation is not big (hard to do in practice without a non-trivial global optimization pass though).

Note that the loop calling log is unrolled twice (hence you should divide the number of instruction by two in the above list to get the number of instructions per call to log).

The target avx2-i32x8 is less impacted by the overheads than the x16 one simply because the avx2-i32x8 code is more latency bound. Thus, the additional instructions can be executed in parallel without adding a significant overhead. Meanwhile, the avx2-i32x16 code already hide latency thanks to double-pumping so additional instruction are not free anymore.

While the generated assembly code can certainly be slightly improved, most of the overheads directly come from the approach/fix so I think there is not much to do (besides using possibly another better approach, if any).

I used my local build but the generated code looks similar to compiler-explorer's code link. I suggest to check which LLVM version you use. Release ISPC and trunk ISPC in compiler-explorer are built with the patched LLVM (see llvm_patches dir in repo). But I don't see difference on my machine between patched LLVM 17 and not-patched LLVM 18 for this example.

Ok. I use the commit 9e74945 @ 20240806 and LLVM 18.1.3. I think I use a non-patched LLVM (from the Ubuntu system packages unless Ubuntu developers patched it).

I measured performance both on GameDev tests and microbenchmarks. For both cases, avx2-i32x8 performance is 20~25 % worse (for tests using log). So, it looks we sacrifice performance for precision for avx2-32x8 also.

Yes. Note that the results are dependent of the target architecture. Out of curiosity, on which architecture/CPU did you run this benchmark?

I measured performance both on GameDev tests and microbenchmarks. For both cases, avx2-i32x8 performance is 20~25 % worse (for tests using log). So, it looks we sacrifice performance for precision for avx2-32x8 also.

Yes. Note that the results are dependent of the target architecture. Out of curiosity, on which architecture/CPU did you run this benchmark?

I tested on Intel i9-12900 (GameDev tests) and AMD 7840HS (microbenchmarks).

I am closing this PR since it is superseded by the recent new one.