[PATCH v1] libm,ieee754:New algorithm of fmod function.

[Wilco Dijkstra]

Hi Kirill,

> I proposed a new algorithm for fmod calculations for double type and
> Algorithm description is in the file.

The new algorithm looks very good and is indeed much faster on 64-bit CPUs!
There are several things that could be done as a followup:

1. Removal of the __fmod veneer which adds an unncessary overhead to set
errno in some cases (this can be done in the special cases similar to how it is
done in eg. exp).

2. In the worst case the new version does ~190 64-bit divisions which has a
significant overhead even with a fast hardware division (I see ~350x difference
between best and worst cases on Neoverse N1). Since all divisions use the same
denominator, it makes sense to use multiplies instead. This would also allow
more bits to be computed per iteration (32 or even 53 if you use a FP divide to
get the 1/y estimate).

> 1) Remove special 64 bit version and use the new version like a default 
> for all platforms.

I'm not convinced by that - 64-bit divisions are slow on 32-bit targets since
they are done in software. A quick check on 32-bit Arm shows it's slower in
cases that need division (the old version is pretty inefficient so could be
improved a lot).

While in general the wordsize-64 version is faster on 32-bit targets, one key
exception is 64-bit division. Using simple iterative division if the exponent
difference is small and multiply by inverse otherwise would be much faster.

> 3) Add benchmark tests for fmod in glibc

In order to make the benchmark useful it shouldn't just have a bunch of random
inputs. Performance varies by more than 2 orders of magnitude across inputs, so
calculating the average over those doesn't give anything useful. The key is to
separately benchmark common cases and worst cases (similar to your github version).
If you look at some of the other inputs, they have several entries like:

## name: slow
## name: subnormal
## name: workload-spec2006.wrf

This allows one to test specific ranges or types of inputs. The "workload-<name>"
variant cycles through the inputs (rather than just repeating the same input many times),
thus stressing branch prediction accurately. The existing workload inputs are extracted
from traces of actual applications, making them even more realistic.

Cheers,
Wilco