New libm fmod function for 64 bit targets.

New libm fmod function for 64 bit targets.

[Kirill Okhotnikov]

Dear Developers,

Thank you very much for the project. I want to contribute to the project 
with a new fmod 64 bit function.

https://sourceware.org/glibc/wiki/Development_Todo/Master#libm_itself

> The C implementations of fmod functions appear suboptimal: they could 
> use __builtin_clz operations for subnormals instead of a loop, where 
> the processor has efficient clz support, and the repeated-subtraction 
> loop may be less efficient than approaches using integer mod (where 
> there is hardware support for that) and repeated squaring, where the 
> exponents are widely different. (Before doing much on this, check what 
> processors would actually benefit.)

During my "research" I found that 32 bit algorithm for fmod written by 
Sun Inc is efficient (*), but it's 64 bit adoption is not. For 64 bit 
system I propose to use a new algorithm which uses integer mod. See 
description here

https://github.com/orex/test_fmod/blob/master/libm-file/e_fmod.c

My tests on x86_64 (Intel, AMD) and ARM64 shows that the new algorithm 
up to 20 times faster for "extreme" cases. And up to two times faster 
for regular using of the function.

https://github.com/orex/test_fmod/blob/master/README.md

Also, I did some unit testing which shows that old and a new algorithms 
gives binary equivalent result for each of billions different pairs (x, 
y) with wide range of numbers including normal, subnormal, and special 
one (NaN INF, 0). Libm tests also passed.

Source code for the tests can be found here.

https://github.com/orex/test_fmod/

The test calculate a hash for all input values and output from libm fmod 
and new fmod. Therefore you can check that it produces the same result 
for different architectures.

Patch is straightforward. You need to change file 
glibc/sysdeps/ieee754/dbl-64/wordsize-64/e_mod.c to 
https://raw.githubusercontent.com/orex/test_fmod/master/libm-file/e_fmod.c.

My question is what else essentially I need to do to apply a patch for 
the library?

Best,
Kirill.

P.S. I develop and wrote the algorithm by myself so the copyright is 
fully mine. Now it is MIT license, but I can change to any license 
needed, of course.

(*) Subnormals loops for CLZ can be changed, of course to builtin_clz. 
But such case looks quite rare in a real calculations.