[PATCH 2/3] math: Optimize dbl-64 remainder implementation

[Wilco Dijkstra]

Hi Adhemerval,

> +    return NAN;
>>
>> Since NAN returns a float NaN, wouldn't C99 nan("") be better?
>
> Right, I think it should ok and not create a nan call for all targets.

We likely need to do something similar as in fmod though, eg. if one is a NaN
return x * y (so invalid operation is set), if y is zero or inf, return __math_edom.

>> Since these branches may be fairly hard to predict, it's worth checking whether
>> this works out better:
>>
>> x = (x >= y) ? x - y : x;
>> x = ((x + x) >= y) ? x - y : x;
>>
>> You can always start with x >= y and only use (x + x) or y * 0.5 for the 2nd part.
>
> I did not see much difference neither on x86_64 (Ryzen 9) or aarch64 (N1) with
> recent gcc.

These are running the same input repeatedly which just tests best-case rather than
average case performance...

I wonder whether we should treat all inputs as workloads and cycle through
the values rather than just repeating?

On 2nd thought, using conditional moves might not get the round to even part of
x - n * y quite right, and there is a risk of triggering inexact in computing x - y early.
I wonder whether we have testcases for these corner cases?

>> In either case, we can remove the if (hy < UINT64_C (0x20000000000000)) by adding
>> an else after the fmod to handle the case of large y.
>
> I am not sure if we can remove it, since it to handle small y.

We can always remove it - small y is only an issue if you use 0.5 * y. If you know y is not
huge (which is true for the call to fmod), then using x + x instead works.

> I think you mean x == -0.0 and yeah I think it should better indeed.

-0.0 and 0.0 are the same - it appears GCC optimizes it to x == 0.0, but it's still wrong to
write x == -0.0 since it gives the incorrect impression that it excludes 0.0 when it does not....

Cheers,
Wilco