[libc][math] fmod/fmodf implementation.

[orex]

Summary

This is a implementation of find remainder fmod function from standard libm.
The underline algorithm is developed by myself, but probably it was first
invented before.
Some features of the implementation:

1. The code is written on more-or-less modern C++.
2. One general implementation for both float and double precision numbers.
3. Spitted platform/architecture dependent and independent code and tests.
4. Tests covers 100% of the code for both float and double numbers. Tests cases with NaN/Inf etc is copied from glibc.
5. The new implementation in general 2-4 times faster for “regular” x,y values. It can be 20 times faster for x/y huge value, but can also be 2 times slower for double denormalized range (according to perf tests provided).
6. Two different implementation of division loop are provided. In some platforms division can be very time consuming operation. Depend on platform it can be 3-10 times slower than multiplication.

Performance tests:

The test is based on core-math project (https://gitlab.inria.fr/core-math/core-math). By Tue Ly suggestion I took hypot function and use it as template for fmod. Preserving all test cases.

./check.sh <--special|--worst> fmodf passed.
CORE_MATH_PERF_MODE=rdtsc ./perf.sh fmodf results are

GNU libc version: 2.35
GNU libc release: stable
21.166 <-- FPU
51.031 <-- current glibc
37.659 <-- this fmod version.

[lntue]

Thanks for adding support for fmod functions with improved performance on regular range! Let's discuss a bit more to see if we can also improve or maintain the performance for denormal inputs.

libc/src/__support/FPUtil/FPBits.h

30

clz are also used in sqrt functions, which would also be used again by FMA function. Would you mind factoring clz functions to another library similar to https://reviews.llvm.org/D124495 ? You can overwrite what I did over there, as this change should be landed before that one. Thanks!

libc/src/__support/FPUtil/generic/CMakeLists.txt

23

Please fix.

libc/src/__support/FPUtil/generic/FMod.h

17

I'm not sure that we can include that many C++ standard headers in here, as it might introduce circular dependency among the libraries. @sivachandra should have known more about these than me.

83

I don't think you would want to use std::isnan or std::isfinite in here. As you can imagine, FPUtil functions would be the ones to provide the *backbone* of libc, and hence std:: functions, so if FPUtil functions depending on std:: math functions, it would likely create circular dependency when building or linking. So it's best to reuse or reimplement those simple math functions in the FPUtil / FPBits themselves.

113

This function seems to better be in FPBits class.

162

This function seems to be generic enough to be a static function in FPBits class.

230

We might not be able to use std::min here.

257

You might have to reimplement std::optional in __support or __support/CPP to prevent circular dependency.

libc/src/math/generic/CMakeLists.txt

1119

Please fix.

libc/test/src/math/differential_testing/CMakeLists.txt

513

Please fix.

libc/test/src/math/exhaustive/FMod_test.cpp

14 ↗

(On Diff #434258)

Is this still needed?

43 ↗

(On Diff #434258)

by * 2^iy is better implemented with ldexp functions https://en.cppreference.com/w/cpp/numeric/math/ldexp

[orex]

Cleanup code by Intue suggestions.

[orex]

Thank you, Intue for your useful comments. I implement all of them.
As for performance issues with denormalized values, I've checked the code of glibc e_fmod.c. It looks like (I assume, but I did not check assembler code) that glibc function comes to denormalaized values faster than in my implementation, where I check most common cases first. I don't think that we should do something with it, because the case is very rare. I can hardly imaging it's real practical usage.

Kirill.

libc/src/__support/FPUtil/FPBits.h

30

I create a new helper with wrapped bulitins. libc/src/__support/FPUtil/builtin_wrappers.h

libc/src/__support/FPUtil/generic/FMod.h

17

All std includes were removed except cerrno.

257

Use "old style" returning optional values. I don't think, that optional reimplementing is needed for this case.

[orex]

Cleanup of the code by Intue suggestions.

[orex]

Cleanup of the code by Intue suggestions

[orex]

Cleanup of the code by Intue suggestions.

[sivachandra]

Thanks a lot for the patch. It seems like it includes a bunch of cleanups / "better this way" items not related to the main goal of the patch. Can you please separate out those parts in to a different patch so that we can keep the review focused? Also, we cannot strictly use std C++ headers. So, no cerrno also. Include errno.h instead.

[orex]

In D127046#3559587, @sivachandra wrote:

Thanks a lot for the patch. It seems like it includes a bunch of cleanups / "better this way" items not related to the main goal of the patch. Can you please separate out those parts in to a different patch so that we can keep the review focused? Also, we cannot strictly use std C++ headers. So, no cerrno also. Include errno.h instead.

Thank you for the comment. I've picked out 3 "better this way" commits:
https://reviews.llvm.org/D127088
https://reviews.llvm.org/D127091
https://reviews.llvm.org/D127097

I'll rebase this commit, when the changes above will be submitted.

[orex]

Rebasing changes on last main.

[orex]

cerrno fix.

[orex]

libc/src/__support/FPUtil/generic/FMod.h

17

cerrno excluded.

83

Implemented using builtin functions.

[orex]

C standard/Posix processing of special numbers.

[lntue]

libc/src/__support/FPUtil/FPBits.h

173

For safety, you might need to add a quick return for when number == 0:

if (unlikely(number == 0)) return zero ...

libc/src/__support/FPUtil/generic/FMod.h

26

nit: treated *separately*

71

There is only one public method for this class, you can replace class with struct and remove public:.

71

Maybe rename the class to FModErrorHandler and add to the comment: following C99 standard with a link to https://en.cppreference.com/w/c/numeric/math/fmod

82

From the C standard, fmod(0, NaN) and fmod(0, inf) should be 0

If x is ±0 and y is not zero, ±0 is returned

88

What's the main purpose of this evaluation?

96

This should be moved to above isnan(x) || isnan(y), I don't think it is hit here.

100

I don't think this line is hit, and all of the unlikely's above are not needed.

104

We don't need this anymore, just use C standard.

145

Instead of using enable_ifwith the boolean InverseMultiplication, it might be better to separate these 2 versions of division_loop into 2 separate structs and then pass them to the second template argument, similar to the way you do exceptional handling. Something like:

// Comments explains what these functions do and what the input parameters are.
template<typename T>
struct FModDivisionLoop {
  static T execute(int n, int hyltzeros, T hx, T hy) { ... }
};
template<typename T>
struct FModInverseMultiplicationLoop {
  static T execute( ... )
};

Then the main class will be something like:

template<typename T, class ErrorHandler = FModErrorHandler<T>, class MainLoop = InverseMultiplicationLoop<T>>
class FMod { ... };

266

FMod<T>::eval() or execute instead of make?

[orex]

libc/src/__support/FPUtil/FPBits.h

173

It was Siva suggestion to move the function here. There is another class called NormalFloat which can handle such things. I propose to move the function back to FMod and implement "full" functionality in that class. What do you think? if number == 0 is not only one case which we need to check. ep value after processing below should also be checked for overflow.

libc/src/__support/FPUtil/generic/FMod.h

26

Sure.

71

Sure.

71

It is not an error handler, but special numbers cases processor. Do you have another name in mind? If no, let's go for Error.

82

Thank you. This is a very good comment. Obviously C standard "F.10.7.1" is not full. It do not describe at all, for example, NaN NaN case. I also do not think, that fmod(0, NaN) should return zero. I prefer to go for better consistent standard which are described here
https://pubs.opengroup.org/onlinepubs/9699919799/functions/fmod.html
https://en.cppreference.com/w/c/numeric/math/fmod
What do you think about this?

88

Thank you for the comment. I don't know any other approach to raise floating-point exception. Please suggest me something better.

96

A question about standard. See above.

100

return false. Yes you are right, but I can't return nothing from the function. It will be a warning warning: non-void function does not return a value in all control paths I don't want my code to produce such things.
unlikely I would like to explicitly say compiler to make just one conditional jump in the function. It is not so big problem, of course, I can remove this.

104

See above.

145

Thank you. Sounds reasonable. I will do this.

266

Yes. Sounds better. Thank you.

[lntue]

libc/src/__support/FPUtil/FPBits.h

173

The main reason you need to check for number == 0 here is because fputil::clz is a wrapper around __built_in_clz and when the inputs are 0, the behavior/output is undefined.

libc/src/__support/FPUtil/generic/FMod.h

71

I think technically these are not errors per se. They actually handle exceptional inputs / outputs, and as a class name, it's more appropriate to use a noun. So I think FModExceptionalHandler or FModExceptionalInputHandler or Helper might be better.

88

You can use https://github.com/llvm/llvm-project/blob/main/libc/src/__support/FPUtil/FEnvImpl.h to set the floating point exceptions. It's a bit more verbose but much clearer. And these are exceptional cases, so we don't care much about the performance anyway.

100

You can remove the last if and put x == 0 in the comment.

[sivachandra]

libc/src/__support/FPUtil/FPBits.h

173

I am unable to locate the context of my suggestion. Can you help me find it?

[orex]

libc/src/__support/FPUtil/FPBits.h

173

Sorry for bothering you. I mixed your suggestion to extract "FBits improvements" and Tue particular suggestion. After new push all old comments are moved to new places for me, so It was difficult for me to figure it out.

173

I never call ctz or clz with zero in fmod case. If the functions now in the separate module, probably we should check zero there? It can have a performance impact, of course, but may be compiler can optimized out. If you always want to check it, I think that it is better to check it there? We can also think about best way to check it. My proposition is to move the function make_value back to FMod and add x == 0 check in ctz/clz functions (probably separate commit).

[lntue]

libc/src/__support/FPUtil/FPBits.h

173

You can use https://reviews.llvm.org/DXXXXXX/new/ instead of https://reviews.llvm.org/DXXXXXX. That should separate old and new comments better.

173

Now that this function make_value is refactored to stay under FPBits class, which might be used beyond fmod. So you will need to either add some comments explicitly stating its non-zero input assumption, or adding an extra check to make sure that it won't silently return wrong outputs when being used by other functions.

[orex]

libc/src/__support/FPUtil/FPBits.h

173

Don't you think, that we need to mention this in clz/ctz builtin wrappers?

[orex]

Some cosmetic changes.

[sivachandra]

libc/src/__support/FPUtil/FPBits.h

173

Don't you think, that we need to mention this in clz/ctz builtin wrappers?

Yes. I think handling the corner case as reasonable should be done within our wrappers. We should not ideally have to pepper user code with checks.

[orex]

libc/src/__support/FPUtil/FPBits.h

173

Thank you for support Siva. I've implemented it in a way it check the case by default.

[lntue]

Thanks for your patient and sorry it took a bit too long for reviewing! Just a few nits and it's good to go.

Also feel free to sync and add fmod to the math status page at libc/docs/math.rst.

libc/src/__support/FPUtil/generic/FMod.h

26

Nit: you don't need to wrap the word separately in *'s, I just use it to highlight my suggestion.

32

Nit: actually without extra conditions on hx, hy, ix, iy, the representation x = hx * 2^ix is ambiguous (i.e., not unique), just from your example. It is only unambiguous if hx, hy, ix, iy are integers and hx, hy are odd, i.e., not divisible by 2.

33–69

So here is my understanding of the algorithm and the math behind it, please correct me if I'm wrong:

The two main properties about the (integer) modulus operation, denoted by mod or %, that we will use are:
For any positive integers a, b, c:

1) a mod b = (a mod (b * c)) mod b
2) (a * c) mod (b * c) = (a mod b) * c

First, let write x = hx * 2^ix and y = hy * 2^iy with hx, hy, ix, iy are integers (assumed to be positive for simplification).
Then the naive implementation of the fmod function with a simple for/while loop:

while (ix > iy) {
  hx = (hx * 2) % hy;
  --ix;
}

is mathematically equivalent to:

x mod y = (hx * 2^ix) mod (hy * 2^iy)
        = ((hx * 2^(ix - iy)) mod hy) * 2^iy     (apply property 2)
        = (( ... (((hx * 2^(ix - iy)) mod (hy * 2^(ix - iy - 1))) mod (hy * 2^(ix - iy - 2)) ... ) mod hy) * 2^iy     (apply property 1 repeatedly)
        = (( ... (( (hx * 2) mod hy) * 2^(ix - iy - 1)) mod (hy * 2^(ix - iy - 2)) ... ) mod hy) * 2^iy     (apply property 2)
        = (( ... (( (hx * 2) mod hy) * 2 ) mod hy ) ... ) * 2) mod hy) * 2^iy     (apply property 2 repeatedly)

And the total number of iterations is ix - iy.

On the other hand, your algorithm exploits the fact that hx, hy are the mantissas of floating point numbers, which use less bits than the storage integers: 24 / 32 for floats and 53 / 64 for doubles, so if in each step of the iteration, we can left shift hx as many bits as the storage integer type can hold, the exponent reduction per step will be at least 32 - 24 = 8 for floats and 64 - 53 = 9 for doubles:

x mod y = (hx * 2^ix) mod (hy * 2^iy)
        = ((hx * 2^(ix - iy) mod hy) * 2^iy
        = (( ... (( (( hx * 2^r1) mod hy ) * 2^r2) mod hy ) ... ) * 2^rk) mod hy) * 2^iy

where r1 + r2 + ... + rk = ix - iy and ri >= sizeof(UInt) - hy length for i = 1..k.

And so the number of iterations is at most by: (ix - iy) / (sizeof(UInt) - mantissa_length).

Feel free to use this to update the comments if it makes sense.

117

exp_diff for n and max_shift for hyltzeroes or other more descriptive names.

179

Maybe using mx, ex or m_x, e_x which are closer to mantissa and exponent of x?

210

Use more descriptive names such as lead_zeros_hy, tail_zeros_hy?

215

max_shift / max_scale_factor maybe more descriptive?

libc/test/src/math/exhaustive/CMakeLists.txt

189

Use lower case for the target name fmod_test.

195

Use lower case for the test file.

[lntue]

libc/src/__support/FPUtil/generic/FMod.h

82

I agree that this makes more sense. We should add a link to https://man7.org/linux/man-pages/man3/fmod.3p.html to support our decision, and maybe send suggestion to cppreference to adjust the order of exceptional cases on their page.

[lntue]

libc/src/__support/FPUtil/generic/FMod.h

237–240

Does marking these 2 conditions unlikely improve throughput?

[lntue]

libc/src/__support/FPUtil/generic/FMod.h

139

Use NumericLimits template from libc/src/__support/CPP/Limits.h.

libc/test/src/math/FModTest.h

37

Use NumericLimits template from libc/src/__support/CPP/Limits.h.

[orex]

Cosmetic changes: veriable renaming, docs update etc.

[orex]

Thank you Tue, for your comments. They was really useful. I'll really appreciate if you go through all my replies. It was a lot of them, so I'm afraid that I can miss something or my changes can be improved even more. I'll also appreciate if you check the performance for fmod in the same way, as other functions. I've attached the file

fmod-core-math.tar.gz52 KBDownload

. You should simply unpack it inside core-math folder.

libc/src/__support/FPUtil/generic/FMod.h

32

Thank you! Just a mistake "unambiguous" -> "ambiguous"

33–69

Thank you! Your explanation is very good. I included some sentences to the description. Huge math inserts looks absolutely unreadable in text format, so I skip them.
Nit: Some preparations step requires for this loop to work. Also this loop can be improved, changing division by subtraction.

117

Done. Thank you!

215

Changed to sides_zeroes_count.

237–240

n == 0 is quite likely condition, I think. But for hx, yes. You are right.

libc/test/src/math/FModTest.h

37

Unfortunately the library is very primitive. It does not have things, I need.

[lntue]

libc/src/__support/FPUtil/FPBits.h

173

Sorry that I confused you and Siva! What I meant was that since the function make_value had undefined behavior when the input is 0, and it was put in the common library, it should either:

• be documented
• add extra checks for 0.

Since you already documented it in item #4, that's good enough for me, and adding extra checks is not needed anymore since there are legitimate usage of make_value and clz when inputs are guaranteed to be non-zero such as in this patch and in other places. You can definitely add extra comments to document the undefined behavior of clz when inputs are 0.

libc/src/__support/FPUtil/builtin_wrappers.h

52

In your use case, extra checks for 0 input are not needed, so commenting on the undefined behavior of builtin_clz when inputs are 0 is good enough. A safe variant with extra checks could be added, but it should use different name. One of the reason is that now if someone one to use builtin_clz without extra checks, they will have to specify the input type T:

clz<T, false>(...)

So when the type T is known, or inside a non-template function, users cannot use type deduction: clz<uint32_t, false>( ... ).
So in summary, my preference is that:

• clz should just be a simple wrap around type matching for __builtin_clz*.
• undefined behavior for inputs 0 can/should be documented.
• safe versions can be added, but use different name so that type deduction can be applied.

libc/test/src/math/FModTest.h

22

You should be able to feed floating point type T directly to expected, and use EXPECT_FQ_EQ macro instead.

37

Look like you only use it for some special constants. With std::numeric_limits, you are actually using <limits> without directly including it which might cause error if some targets or configs that do not transitively include it.

It would be better to add those constants that you need to constant creating functions in FPBits class, and/or update the DECLARE_SPECIAL_CONSTANTS macro at https://github.com/llvm/llvm-project/blob/main/libc/utils/UnitTest/FPMatcher.h#L70 . It could a separate patch that this one can depend on.

[orex]

Tests polishing.

[orex]

Tue. I think, if we will not find agreement with ctz/clz syntax, we can ask Siva to solve the problem. Also, as long as you will use the function in the next changes, you can also try them and if it will not work good, improve it.

libc/src/__support/FPUtil/builtin_wrappers.h

52

I have a little bit different opinion for this case. If somebody would like to simply use clz/ctz he can use it without any complication and in a safe way. But, if somebody know, that he can improve the performance, skipping the test, it looks like OK to have such complications. For example. I don't see much difference between, adding new function have syntax like this ctz<decltype(i), false>(i). Moreover such syntax can help to avoid implicit type conversion. Another point to use such template is to pass check zero parameter above easily. For example, make_value function can looks like make_value<bool CheckZero>(.... Of course it can be done with if consexpr, but template will be more straightforward.

libc/test/src/math/FModTest.h

22

Thank you.

37

Can you explain your point, please. From my point of view, I just forget to #include <limits>. From another side, I have a feeling, from you comment, that I should not do this. Is it so? If yes, why. It is a test. Tests are "final" instances, so they can include whatever they want. Fre example, I've checked exhaustive_test.cpp. It includes "half of" STL.

[lntue]

libc/src/__support/FPUtil/builtin_wrappers.h

52

So from the readability standpoint of both implementers and reviewers, I think clz(x) and safe_clz(x) or clz(x) and unsafe_clz(x)' are better to understand than clz(x) and clz<delctype(x), false>(x)`.
I'm also not a big fan of using boolean in template parameter, unless the template name makes it very clear what true and false mean. If it's only used in one or a few tests, that's fine. But this is on a utility header that will be used everywhere. If we really want to pass these as template parameters, it's better to wrap them in functors, like function<SafeClz> vs function<true>. Hope it makes sense?

Also, by changing the default clz(x) in this patch, you would need to update all current usages of it like in sqrt to the default version.

libc/test/src/math/FModTest.h

37

For unit tests, we try to limits the use of std:: and C++ standard header. You can see other unit tests in libc/test/src/math. Exhaustive tests are a bit different. They are kind of integration tests, so we do not be so strict about that (yet). So for example, std::numeric_limit<>::quiet_NaN() (most likely) simply call __builtin_nan*, which technically what our library implements explicitly. You don't have to do it now because this test currently does not have that problem. But we should clean it up in the future if not now.

[orex]

unsafe ctz/clz

[orex]

libc/src/__support/FPUtil/builtin_wrappers.h

52

OK. You won.)))) I've change all occurrence to unsafe_ctz/clz.

[lntue]

Thanks for sticking with me until now! Please sync to head and feel free to land when the pre-merge checks turn green.

[orex]

Changed unsafe_clz to safe_clz in string_to_float.h

libc/config/darwin/arm/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.fma
     libc.src.math.fmaf
     libc.src.math.fmax
     libc.src.math.fmaxf
     libc.src.math.fmaxl
     libc.src.math.fmin
     libc.src.math.fminf
     libc.src.math.fminl
+    libc.src.math.fmod
+    libc.src.math.fmodf
     libc.src.math.frexp
     libc.src.math.frexpf
     libc.src.math.frexpl
     libc.src.math.hypot
     libc.src.math.hypotf
     libc.src.math.ilogb
     libc.src.math.ilogbf
     libc.src.math.ilogbl

libc/config/linux/aarch64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.fma
     libc.src.math.fmaf
     libc.src.math.fmax
     libc.src.math.fmaxf
     libc.src.math.fmaxl
     libc.src.math.fmin
     libc.src.math.fminf
     libc.src.math.fminl
+    libc.src.math.fmod
+    libc.src.math.fmodf
     libc.src.math.frexp
     libc.src.math.frexpf
     libc.src.math.frexpl
     libc.src.math.hypot
     libc.src.math.hypotf
     libc.src.math.ilogb
     libc.src.math.ilogbf
     libc.src.math.ilogbl

libc/config/linux/x86_64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.fma
     libc.src.math.fmaf
     libc.src.math.fmin
     libc.src.math.fminf
     libc.src.math.fminl
     libc.src.math.fmax
     libc.src.math.fmaxf
     libc.src.math.fmaxl
+    libc.src.math.fmod
+    libc.src.math.fmodf
     libc.src.math.frexp
     libc.src.math.frexpf
     libc.src.math.frexpl
     libc.src.math.hypot
     libc.src.math.hypotf
     libc.src.math.ilogb
     libc.src.math.ilogbf
     libc.src.math.ilogbl

libc/config/windows/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.fma
     libc.src.math.fmaf
     libc.src.math.fmin
     libc.src.math.fminf
     libc.src.math.fminl
     libc.src.math.fmax
     libc.src.math.fmaxf
     libc.src.math.fmaxl
+    libc.src.math.fmod
+    libc.src.math.fmodf
     libc.src.math.frexp
     libc.src.math.frexpf
     libc.src.math.frexpl
     libc.src.math.hypot
     libc.src.math.hypotf
     libc.src.math.ilogb
     libc.src.math.ilogbf
     libc.src.math.ilogbl

libc/docs/math.rst

 ============== ================ =============== ======================
 ceil           |check|          |check|         |check|  
 copysign       |check|          |check|         |check|  
 fabs           |check|          |check|         |check|  
 fdim           |check|          |check|         |check|  
 floor          |check|          |check|         |check|  
 fmax           |check|          |check|         |check|  
 fmin           |check|          |check|         |check|  
-fmod
+fmod           |check|          |check|
 fpclassify
 frexp          |check|          |check|         |check|  
 ilogb          |check|          |check|         |check|  
 isfinite
 isgreater
 isgreaterequal
 isinf
 isless

libc/spec/stdc.td

HeaderSpec Math = HeaderSpec<
 
           FunctionSpec<"fmax", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
           FunctionSpec<"fmaxf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
           FunctionSpec<"fmaxl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoubleType>]>,
 
           FunctionSpec<"fma", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
           FunctionSpec<"fmaf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>, ArgSpec<FloatType>]>,
 
+          FunctionSpec<"fmod", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
+
+          FunctionSpec<"fmodf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
+
           FunctionSpec<"frexp", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<IntPtr>]>,
           FunctionSpec<"frexpf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<IntPtr>]>,
           FunctionSpec<"frexpl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<IntPtr>]>,
 
           FunctionSpec<"hypot", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
           FunctionSpec<"hypotf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
 
           FunctionSpec<"ilogb", RetValSpec<IntType>, [ArgSpec<DoubleType>]>,

libc/src/__support/FPUtil/FPBits.h

 //===-- Abstract class for bit manipulation of float numbers. ---*- C++ -*-===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #ifndef LLVM_LIBC_SRC_SUPPORT_FPUTIL_FP_BITS_H
 #define LLVM_LIBC_SRC_SUPPORT_FPUTIL_FP_BITS_H
 
 #include "PlatformDefs.h"
 
 #include "src/__support/CPP/Bit.h"
 #include "src/__support/CPP/TypeTraits.h"
+#include "src/__support/FPUtil/builtin_wrappers.h"
+#include "src/__support/common.h"
 
 #include "FloatProperties.h"
 #include <stdint.h>
 
 namespace __llvm_libc {
 namespace fputil {
 
 template <typename T> struct MantissaWidth {
   static constexpr unsigned VALUE = FloatProperties<T>::MANTISSA_WIDTH;
 };
 
 template <typename T> struct ExponentWidth {
   static constexpr unsigned VALUE = FloatProperties<T>::EXPONENT_WIDTH;
 };
 
 // A generic class to represent single precision, double precision, and quad
 // precision IEEE 754 floating point formats.
 // On most platforms, the 'float' type corresponds to single precision floating
 // point numbers, the 'double' type corresponds to double precision floating
 // point numers, and the 'long double' type corresponds to the quad precision
 // floating numbers. On x86 platforms however, the 'long double' type maps to
 // an x87 floating point format. This format is an IEEE 754 extension format.
 // It is handled as an explicit specialization of this class.
static constexpr FPBits<T> neg_inf() {
     return bits;
   }
 
   static constexpr T build_nan(UIntType v) {
     FPBits<T> bits = inf();
     bits.set_mantissa(v);
     return T(bits);
   }
+
+  // The function convert integer number and unbiased exponent to proper float
+  // T type:
+  //   Result = number * 2^(ep+1 - exponent_bias)
+  // Be careful!
+  //   1) "ep" is raw exponent value.
+  //   2) The function add to +1 to ep for seamless normalized to denormalized
+  //      transition.
+  //   3) The function did not check exponent high limit.
+  //   4) "number" zero value is not processed correctly.
+  //   5) Number is unsigned, so the result can be only positive.
+  inline static constexpr FPBits<T> make_value(UIntType number, int ep) {
+    FPBits<T> result;
+    // offset: +1 for sign, but -1 for implicit first bit
+    int lz = fputil::unsafe_clz(number) - FloatProp::EXPONENT_WIDTH;
+    number <<= lz;
+    ep -= lz;
+
+    if (likely(ep >= 0)) {
+      // Implicit number bit will be removed by mask
+      result.set_mantissa(number);
+      result.set_unbiased_exponent(ep + 1);
+    } else {
+      result.set_mantissa(number >> -ep);
+    }
+    return result;
+  }
 };
 
 } // namespace fputil
 } // namespace __llvm_libc
 
 #ifdef SPECIAL_X86_LONG_DOUBLE
 #include "x86_64/LongDoubleBits.h"
 #endif
 
 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_FP_BITS_H

libc/src/__support/FPUtil/Hypot.h

 namespace fputil {
 
 namespace internal {
 
 template <typename T>
 static inline T find_leading_one(T mant, int &shift_length) {
   shift_length = 0;
   if (mant > 0) {
-    shift_length = (sizeof(mant) * 8) - 1 - clz(mant);
+    shift_length = (sizeof(mant) * 8) - 1 - unsafe_clz(mant);
   }
   return T(1) << shift_length;
 }
 
 } // namespace internal
 
 template <typename T> struct DoubleLength;
 

libc/src/__support/FPUtil/builtin_wrappers.h

 
 namespace __llvm_libc {
 namespace fputil {
 
 // The following overloads are matched based on what is accepted by
 // __builtin_clz/ctz* rather than using the exactly-sized aliases from stdint.h.
 // This way, we can avoid making any assumptions about integer sizes and let the
 // compiler match for us.
+namespace __internal {
+
+template <typename T> static inline int correct_zero(T val, int bits) {
+  if (val == T(0))
+    return sizeof(T(0)) * 8;
+  else
+    return bits;
+}
+
 template <typename T> static inline int clz(T val);
 template <> inline int clz<unsigned int>(unsigned int val) {
   return __builtin_clz(val);
 }
 template <> inline int clz<unsigned long int>(unsigned long int val) {
   return __builtin_clzl(val);
 }
 template <> inline int clz<unsigned long long int>(unsigned long long int val) {
   return __builtin_clzll(val);
 }
 
 template <typename T> static inline int ctz(T val);
 template <> inline int ctz<unsigned int>(unsigned int val) {
   return __builtin_ctz(val);
 }
 template <> inline int ctz<unsigned long int>(unsigned long int val) {
   return __builtin_ctzl(val);
 }
 template <> inline int ctz<unsigned long long int>(unsigned long long int val) {
   return __builtin_ctzll(val);
 }
+} // namespace __internal
+
+template <typename T> static inline int safe_ctz(T val) {
+  return __internal::correct_zero(val, __internal::ctz(val));
+}
+
+template <typename T> static inline int unsafe_ctz(T val) {
+  return __internal::ctz(val);
+}
+
+template <typename T> static inline int safe_clz(T val) {
+  return __internal::correct_zero(val, __internal::clz(val));
+}
+
+template <typename T> static inline int unsafe_clz(T val) {
+  return __internal::clz(val);
+}
 
 template <typename T> static inline bool isnan(T val) {
   return __builtin_isnan(val);
 }
 
 template <typename T> static inline bool isinf(T val) {
   return __builtin_isinf(val);
 }

libc/src/__support/FPUtil/generic/CMakeLists.txt

 add_header_library(
   sqrt
   HDRS
     sqrt.h
     sqrt_80_bit_long_double.h
   DEPENDS
     libc.src.__support.CPP.uint128
 )
 
 add_header_library(
   fma
   HDRS
     FMA.h
   DEPENDS
     libc.src.__support.CPP.uint128
 )
+
+add_header_library(
+  fmod
+  HDRS
+    FMod.h
+)

libc/src/__support/FPUtil/generic/FMA.h

template <> inline double fma<double>(double x, double y, double z) {
   }
 
   uint64_t result = 0;
   int r_exp = 0; // Unbiased exponent of the result
 
   // Normalize the result.
   if (prod_mant != 0) {
     uint64_t prod_hi = static_cast<uint64_t>(prod_mant >> 64);
-    int lead_zeros =
-        prod_hi ? clz(prod_hi) : 64 + clz(static_cast<uint64_t>(prod_mant));
+    int lead_zeros = prod_hi
+                         ? unsafe_clz(prod_hi)
+                         : 64 + unsafe_clz(static_cast<uint64_t>(prod_mant));
     // Move the leading 1 to the most significant bit.
     prod_mant <<= lead_zeros;
     // The lower 64 bits are always sticky bits after moving the leading 1 to
     // the most significant bit.
     sticky_bits |= (static_cast<uint64_t>(prod_mant) != 0);
     result = static_cast<uint64_t>(prod_mant >> 64);
     // Change prod_lsb_exp the be the exponent of the least significant bit of
     // the result.

libc/src/__support/FPUtil/generic/FMod.h

+//===-- Common header for fmod implementations ------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_SUPPORT_FPUTIL_GENERIC_FMOD_H
+#define LLVM_LIBC_SRC_SUPPORT_FPUTIL_GENERIC_FMOD_H
+
+#include "src/__support/CPP/Limits.h"
+#include "src/__support/CPP/TypeTraits.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/builtin_wrappers.h"
+#include "src/__support/common.h"
+#include "src/math/generic/math_utils.h"
+
+namespace __llvm_libc {
+namespace fputil {
+namespace generic {
+
+//  Objective:
+//    The  algorithm uses  integer arithmetic  (max uint64_t)  for general case.
+//    Some common  cases, like  abs(x) < abs(y)  or  abs(x) < 1000 *  abs(y) are
+//    treated specially to increase  performance.  The part of checking  special
+//    cases, numbers NaN, INF etc. treated separately.
+//
+//  Objective:
+//    1) FMod definition (https://cplusplus.com/reference/cmath/fmod/):
+//       fmod = numer - tquot * denom, where tquot is the truncated
+//       (i.e., rounded towards zero) result of: numer/denom.
+//    2) FMod with negative x and/or y can be trivially converted to fmod for
+//       positive x and y. Therefore the algorithm below works only with
+//       positive numbers.
+//    3) All positive floating point numbers can be represented as m * 2^e,
+//       where "m" is positive integer and "e" is signed.
+//    4) FMod function can be calculated in integer numbers (x > y):
+//         fmod = m_x * 2^e_x - tquot * m_y * 2^e_y
+//              = 2^e_y * (m_x * 2^(e_x - e^y) - tquot * m_y).
+//       All variables in parentheses are unsigned integers.
+//
+//  Mathematical background:
+//    Input x,y in the algorithm is represented (mathematically) like m_x*2^e_x
+//    and m_y*2^e_y. This is an ambiguous number representation. For example:
+//      m * 2^e = (2 * m) * 2^(e-1)
+//    The algorithm uses the facts that
+//      r = a % b = (a % (N * b)) % b,
+//      (a * c) % (b * c) = (a % b) * c
+//    where N is positive  integer number. a, b and c - positive. Let's  adopt
+//    the formula for representation above.
+//      a = m_x * 2^e_x, b = m_y * 2^e_y, N = 2^k
+//      r(k) = a % b = (m_x * 2^e_x) % (2^k * m_y * 2^e_y)
+//           = 2^(e_y + k) * (m_x * 2^(e_x - e_y - k) % m_y)
+//      r(k) = m_r * 2^e_r = (m_x % m_y) * 2^(m_y + k)
+//           = (2^p * (m_x % m_y) * 2^(e_y + k - p))
+//        m_r = 2^p * (m_x % m_y), e_r = m_y + k - p
+//
+//  Algorithm description:
+//  First, let write x = m_x * 2^e_x and y = m_y * 2^e_y with m_x, m_y, e_x, e_y
+//  are integers (m_x amd m_y positive).
+//  Then the naive  implementation of the fmod function with a simple
+//  for/while loop:
+//      while (e_x > e_y) {
+//        m_x *= 2; --e_x; //  m_x * 2^e_x == 2 * m_x * 2^(e_x - 1)
+//        m_x %= m_y;
+//      }
+//  On the other hand, the algorithm exploits the fact that m_x, m_y are the
+//  mantissas of floating point numbers, which use less bits than the storage
+//  integers: 24 / 32 for floats and 53 / 64 for doubles, so if in each step of
+//  the iteration, we can left shift m_x as many bits as the storage integer
+//  type can hold, the exponent reduction per step will be at least 32 - 24 = 8
+//  for floats and 64 - 53 = 11 for doubles (double example below):
+//      while (e_x > e_y) {
+//        m_x <<= 11; e_x -= 11; //  m_x * 2^e_x == 2^11 * m_x * 2^(e_x - 11)
+//        m_x %= m_y;
+//      }
+//  Some extra improvements are done:
+//    1) Shift m_y maximum to the right, which can significantly improve
+//       performance for small integer numbers (y = 3 for example).
+//       The m_x shift in the loop can be 62 instead of 11 for double.
+//    2) For some architectures with very slow division, it can be better to
+//       calculate inverse value ones, and after do multiplication in the loop.
+//    3) "likely" special cases are treated specially to improve performance.
+//
+//  Simple example:
+//  The examples below use byte for simplicity.
+//    1) Shift hy maximum to right without losing bits and increase iy value
+//       m_y = 0b00101100 e_y = 20 after shift m_y = 0b00001011 e_y = 22.
+//    2) m_x = m_x % m_y.
+//    3) Move m_x maximum to left. Note that after (m_x = m_x % m_y) CLZ in m_x
+//    is not lower than CLZ in m_y. m_x=0b00001001 e_x = 100, m_x=0b10010000,
+//       e_x = 100-4 = 96.
+//    4) Repeat (2) until e_x == e_y.
+//
+//  Complexity analysis (double):
+//    Converting x,y to (m_x,e_x),(m_y, e_y): CTZ/shift/AND/OR/if. Loop  count:
+//      (m_x - m_y) / (64 -  "length of m_y").
+//      max("length of m_y")  = 53,
+//      max(e_x - e_y)  = 2048
+//    Maximum operation is  186. For rare "unrealistic" cases.
+//
+//  Special cases (double):
+//    Supposing  that  case  where |y| > 1e-292 and |x/y|<2000  is  very  common
+//    special processing is implemented. No m_y alignment, no loop:
+//      result = (m_x * 2^(e_x - e_y)) % m_y.
+//    When x and y are both subnormal (rare case but...) the
+//      result = m_x % m_y.
+//    Simplified conversion back to double.
+
+// Exceptional cases handler according to cppreference.com
+//    https://en.cppreference.com/w/cpp/numeric/math/fmod
+// and POSIX standard described in Linux man
+//   https://man7.org/linux/man-pages/man3/fmod.3p.html
+// C standard for the function is not full, so not by default (although it can
+// be implemented in another handler.
+template <typename T> struct FModExceptionalInputHandler {
+
+  static_assert(cpp::IsFloatingPointType<T>::Value,
+                "FModCStandardWrapper instantiated with invalid type.");
+
+  static bool PreCheck(T x, T y, T &out) {
+    if (likely(y != 0 && fputil::isfinite(y) && fputil::isfinite(x))) {
+      return false;
+    }
+
+    if (fputil::isnan(x) || fputil::isnan(y)) {
+      out = fputil::quiet_NaN(T(0));
+      return true;
+    }
+
+    if (fputil::isinf(x) || y == 0) {
+      fputil::set_except(FE_INVALID);
+      out = with_errno(fputil::quiet_NaN(T(0)), EDOM);
+      return true;
+    }
+
+    if (fputil::isinf(y)) {
+      out = x;
+      return true;
+    }
+
+    // case where x == 0
+    out = x;
+    return true;
+  }
+};
+
+template <typename T> struct FModFastMathWrapper {
+
+  static_assert(cpp::IsFloatingPointType<T>::Value,
+                "FModFastMathWrapper instantiated with invalid type.");
+
+  static bool PreCheck(T, T, T &) { return false; }
+};
+
+template <typename T> class FModDivisionSimpleHelper {
+private:
+  using intU_t = typename FPBits<T>::UIntType;
+
+public:
+  inline constexpr static intU_t execute(int exp_diff, int sides_zeroes_count,
+                                         intU_t m_x, intU_t m_y) {
+    while (exp_diff > sides_zeroes_count) {
+      exp_diff -= sides_zeroes_count;
+      m_x <<= sides_zeroes_count;
+      m_x %= m_y;
+    }
+    m_x <<= exp_diff;
+    m_x %= m_y;
+    return m_x;
+  }
+};
+
+template <typename T> class FModDivisionInvMultHelper {
+private:
+  using FPB = FPBits<T>;
+  using intU_t = typename FPB::UIntType;
+
+public:
+  inline constexpr static intU_t execute(int exp_diff, int sides_zeroes_count,
+                                         intU_t m_x, intU_t m_y) {
+    if (exp_diff > sides_zeroes_count) {
+      intU_t inv_hy = (cpp::NumericLimits<intU_t>::max() / m_y);
+      while (exp_diff > sides_zeroes_count) {
+        exp_diff -= sides_zeroes_count;
+        intU_t hd =
+            (m_x * inv_hy) >> (FPB::FloatProp::BIT_WIDTH - sides_zeroes_count);
+        m_x <<= sides_zeroes_count;
+        m_x -= hd * m_y;
+        while (unlikely(m_x > m_y))
+          m_x -= m_y;
+      }
+      intU_t hd = (m_x * inv_hy) >> (FPB::FloatProp::BIT_WIDTH - exp_diff);
+      m_x <<= exp_diff;
+      m_x -= hd * m_y;
+      while (unlikely(m_x > m_y))
+        m_x -= m_y;
+    } else {
+      m_x <<= exp_diff;
+      m_x %= m_y;
+    }
+    return m_x;
+  }
+};
+
+template <typename T, class Wrapper = FModExceptionalInputHandler<T>,
+          class DivisionHelper = FModDivisionSimpleHelper<T>>
+class FMod {
+  static_assert(cpp::IsFloatingPointType<T>::Value,
+                "FMod instantiated with invalid type.");
+
+private:
+  using FPB = FPBits<T>;
+  using intU_t = typename FPB::UIntType;
+
+  inline static constexpr FPB eval_internal(FPB sx, FPB sy) {
+
+    if (likely(sx.uintval() <= sy.uintval())) {
+      if (sx.uintval() < sy.uintval())
+        return sx;        // |x|<|y| return x
+      return FPB::zero(); // |x|=|y| return 0.0
+    }
+
+    int e_x = sx.get_unbiased_exponent();
+    int e_y = sy.get_unbiased_exponent();
+
+    // Most common case where |y| is "very normal" and |x/y| < 2^EXPONENT_WIDTH
+    if (likely(e_y > int(FPB::FloatProp::MANTISSA_WIDTH) &&
+               e_x - e_y <= int(FPB::FloatProp::EXPONENT_WIDTH))) {
+      intU_t m_x = sx.get_explicit_mantissa();
+      intU_t m_y = sy.get_explicit_mantissa();
+      intU_t d = (e_x == e_y) ? (m_x - m_y) : (m_x << (e_x - e_y)) % m_y;
+      if (d == 0)
+        return FPB::zero();
+      // iy - 1 because of "zero power" for number with power 1
+      return FPB::make_value(d, e_y - 1);
+    }
+    /* Both subnormal special case. */
+    if (unlikely(e_x == 0 && e_y == 0)) {
+      FPB d;
+      d.set_mantissa(sx.uintval() % sy.uintval());
+      return d;
+    }
+
+    // Note that hx is not subnormal by conditions above.
+    intU_t m_x = sx.get_explicit_mantissa();
+    e_x--;
+
+    intU_t m_y = sy.get_explicit_mantissa();
+    int lead_zeros_m_y = FPB::FloatProp::EXPONENT_WIDTH;
+    if (likely(e_y > 0)) {
+      e_y--;
+    } else {
+      m_y = sy.get_mantissa();
+      lead_zeros_m_y = unsafe_clz(m_y);
+    }
+
+    // Assume hy != 0
+    int tail_zeros_m_y = unsafe_ctz(m_y);
+    int sides_zeroes_count = lead_zeros_m_y + tail_zeros_m_y;
+    // n > 0 by conditions above
+    int exp_diff = e_x - e_y;
+    {
+      // Shift hy right until the end or n = 0
+      int right_shift = exp_diff < tail_zeros_m_y ? exp_diff : tail_zeros_m_y;
+      m_y >>= right_shift;
+      exp_diff -= right_shift;
+      e_y += right_shift;
+    }
+
+    {
+      // Shift hx left until the end or n = 0
+      int left_shift = exp_diff < int(FPB::FloatProp::EXPONENT_WIDTH)
+                           ? exp_diff
+                           : FPB::FloatProp::EXPONENT_WIDTH;
+      m_x <<= left_shift;
+      exp_diff -= left_shift;
+    }
+
+    m_x %= m_y;
+    if (unlikely(m_x == 0))
+      return FPB::zero();
+
+    if (exp_diff == 0)
+      return FPB::make_value(m_x, e_y);
+
+    /* hx next can't be 0, because hx < hy, hy % 2 == 1 hx * 2^i % hy != 0 */
+    m_x = DivisionHelper::execute(exp_diff, sides_zeroes_count, m_x, m_y);
+    return FPB::make_value(m_x, e_y);
+  }
+
+public:
+  static inline T eval(T x, T y) {
+    if (T out; Wrapper::PreCheck(x, y, out))
+      return out;
+    FPB sx(x), sy(y);
+    bool sign = sx.get_sign();
+    sx.set_sign(false);
+    sy.set_sign(false);
+    FPB result = eval_internal(sx, sy);
+    result.set_sign(sign);
+    return result.get_val();
+  }
+};
+
+} // namespace generic
+} // namespace fputil
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_GENERIC_FMOD_H

libc/src/__support/FPUtil/generic/sqrt.h

 template <> struct SpecialLongDouble<long double> {
   static constexpr bool VALUE = true;
 };
 #endif // SPECIAL_X86_LONG_DOUBLE
 
 template <typename T>
 static inline void normalize(int &exponent,
                              typename FPBits<T>::UIntType &mantissa) {
-  const int shift =
-      clz(mantissa) - (8 * sizeof(mantissa) - 1 - MantissaWidth<T>::VALUE);
+  const int shift = unsafe_clz(mantissa) -
+                    (8 * sizeof(mantissa) - 1 - MantissaWidth<T>::VALUE);
   exponent -= shift;
   mantissa <<= shift;
 }
 
 #ifdef LONG_DOUBLE_IS_DOUBLE
 template <>
 inline void normalize<long double>(int &exponent, uint64_t &mantissa) {
   normalize<double>(exponent, mantissa);

libc/src/__support/FPUtil/generic/sqrt_80_bit_long_double.h

 #include "src/__support/FPUtil/builtin_wrappers.h"
 
 namespace __llvm_libc {
 namespace fputil {
 namespace x86 {
 
 inline void normalize(int &exponent, UInt128 &mantissa) {
   const int shift =
-      clz(static_cast<uint64_t>(mantissa)) -
+      unsafe_clz(static_cast<uint64_t>(mantissa)) -
       (8 * sizeof(uint64_t) - 1 - MantissaWidth<long double>::VALUE);
   exponent -= shift;
   mantissa <<= shift;
 }
 
 // if constexpr statement in sqrt.h still requires x86::sqrt to be declared
 // even when it's not used.
 static inline long double sqrt(long double x);

libc/src/__support/str_to_float.h

template <class T> uint32_t inline leading_zeroes(T inputNumber) {
   }
   if (inputNumber >> cur_guess > 0) {
     cur_guess++;
   }
   return BITS_IN_T - cur_guess;
 }
 
 template <> uint32_t inline leading_zeroes<uint32_t>(uint32_t inputNumber) {
-  return inputNumber == 0 ? 32 : fputil::clz(inputNumber);
+  return fputil::safe_clz(inputNumber);
 }
 
 template <> uint32_t inline leading_zeroes<uint64_t>(uint64_t inputNumber) {
-  return inputNumber == 0 ? 64 : fputil::clz(inputNumber);
+  return fputil::safe_clz(inputNumber);
 }
 
 static inline uint64_t low64(const UInt128 &num) {
   return static_cast<uint64_t>(num & 0xffffffffffffffff);
 }
 
 static inline uint64_t high64(const UInt128 &num) {
   return static_cast<uint64_t>(num >> 64);

libc/src/math/CMakeLists.txt

 add_math_entrypoint_object(fmax)
 add_math_entrypoint_object(fmaxf)
 add_math_entrypoint_object(fmaxl)
 
 add_math_entrypoint_object(fmin)
 add_math_entrypoint_object(fminf)
 add_math_entrypoint_object(fminl)
 
+add_math_entrypoint_object(fmod)
+add_math_entrypoint_object(fmodf)
+
 add_math_entrypoint_object(frexp)
 add_math_entrypoint_object(frexpf)
 add_math_entrypoint_object(frexpl)
 
 add_math_entrypoint_object(hypot)
 add_math_entrypoint_object(hypotf)
 
 add_math_entrypoint_object(ilogb)

libc/src/math/fmod.h

+//===-- Implementation header for fmod --------------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_FMOD_H
+#define LLVM_LIBC_SRC_MATH_FMOD_H
+
+namespace __llvm_libc {
+
+double fmod(double x, double y);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_FMOD_H

libc/src/math/fmodf.h

+//===-- Implementation header for fmodf -------------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_FMODF_H
+#define LLVM_LIBC_SRC_MATH_FMODF_H
+
+namespace __llvm_libc {
+
+float fmodf(float x, float y);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_FMODF_H

libc/src/math/generic/CMakeLists.txt

HDRS
     dp_trig.h
   DEPENDS
     libc.src.__support.FPUtil.fputil #FPBits and ManipulationFunction
     libc.src.__support.FPUtil.xfloat
     libc.src.__support.CPP.uint
   COMPILE_OPTIONS
    -O3
 )
+
+add_entrypoint_object(
+  fmod
+  SRCS
+    fmod.cpp
+  HDRS
+    ../fmod.h
+  DEPENDS
+    libc.include.math
+    libc.src.__support.FPUtil.generic.fmod
+  COMPILE_OPTIONS
+    -O3
+)
+
+add_entrypoint_object(
+  fmodf
+  SRCS
+    fmodf.cpp
+  HDRS
+    ../fmodf.h
+  DEPENDS
+    libc.include.math
+    libc.src.__support.FPUtil.generic.fmod
+  COMPILE_OPTIONS
+    -O3
+)
 No newline at end of file

libc/src/math/generic/fmod.cpp

+//===-- Double-precision fmod function ------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/fmod.h"
+#include "src/__support/FPUtil/generic/FMod.h"
+#include "src/__support/common.h"
+
+namespace __llvm_libc {
+
+LLVM_LIBC_FUNCTION(double, fmod, (double x, double y)) {
+  return fputil::generic::FMod<double>::eval(x, y);
+}
+
+} // namespace __llvm_libc

libc/src/math/generic/fmodf.cpp

+//===-- Single-precision fmodf function -----------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/fmodf.h"
+#include "src/__support/FPUtil/generic/FMod.h"
+#include "src/__support/common.h"
+
+namespace __llvm_libc {
+
+LLVM_LIBC_FUNCTION(float, fmodf, (float x, float y)) {
+  return fputil::generic::FMod<float>::eval(x, y);
+}
+
+} // namespace __llvm_libc

libc/test/src/math/CMakeLists.txt

add_fp_unittest(
   SRCS
     log1pf_test.cpp
   DEPENDS
     libc.include.math
     libc.src.math.log1pf
     libc.src.__support.FPUtil.fputil
 )
 
+add_fp_unittest(
+  fmodf_test
+  SUITE
+    libc_math_unittests
+  SRCS
+    fmodf_test.cpp
+  HDRS
+    FModTest.h
+  DEPENDS
+    libc.include.math
+    libc.src.math.fmodf
+    libc.src.__support.FPUtil.fputil
+)
+
+add_fp_unittest(
+  fmod_test
+  SUITE
+    libc_math_unittests
+  SRCS
+    fmod_test.cpp
+  HDRS
+    FModTest.h
+  DEPENDS
+    libc.include.math
+    libc.src.math.fmod
+    libc.src.__support.FPUtil.fputil
+)
+
 add_subdirectory(generic)
 add_subdirectory(exhaustive)
 add_subdirectory(differential_testing)

libc/test/src/math/FModTest.h

+//===-- Utility class to test fmod special numbers ------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_TEST_SRC_MATH_FMODTEST_H
+#define LLVM_LIBC_TEST_SRC_MATH_FMODTEST_H
+
+#include "src/__support/FPUtil/BasicOperations.h"
+#include "src/__support/FPUtil/NearestIntegerOperations.h"
+#include "utils/UnitTest/FPMatcher.h"
+#include "utils/UnitTest/Test.h"
+
+#include <limits>
+#include <math.h>
+
+#define TEST_SPECIAL(x, y, expected, dom_err, expected_exception)              \
+  EXPECT_FP_EQ(expected, f(x, y));                                             \
+  EXPECT_MATH_ERRNO((dom_err) ? EDOM : 0);                                     \
+  EXPECT_FP_EXCEPTION(expected_exception);                                     \
+  __llvm_libc::fputil::clear_except(FE_ALL_EXCEPT)
+
+#define TEST_REGULAR(x, y, expected) TEST_SPECIAL(x, y, expected, false, 0)
+
+template <typename T> class FmodTest : public __llvm_libc::testing::Test {
+
+  DECLARE_SPECIAL_CONSTANTS(T)
+
+public:
+  typedef T (*FModFunc)(T, T);
+
+  void testSpecialNumbers(FModFunc f) {
+    using nl = std::numeric_limits<T>;
+
+    // fmod (+0, y) == +0 for y != 0.
+    TEST_SPECIAL(0.0, 3.0, 0.0, false, 0);
+    TEST_SPECIAL(0.0, nl::denorm_min(), 0.0, false, 0);
+    TEST_SPECIAL(0.0, -nl::denorm_min(), 0.0, false, 0);
+    TEST_SPECIAL(0.0, nl::min(), 0.0, false, 0);
+    TEST_SPECIAL(0.0, -nl::min(), 0.0, false, 0);
+    TEST_SPECIAL(0.0, nl::max(), 0.0, false, 0);
+    TEST_SPECIAL(0.0, -nl::max(), 0.0, false, 0);
+
+    // fmod (-0, y) == -0 for y != 0.
+    TEST_SPECIAL(neg_zero, 3.0, neg_zero, false, 0);
+    TEST_SPECIAL(neg_zero, nl::denorm_min(), neg_zero, false, 0);
+    TEST_SPECIAL(neg_zero, -nl::denorm_min(), neg_zero, false, 0);
+    TEST_SPECIAL(neg_zero, nl::min(), neg_zero, false, 0);
+    TEST_SPECIAL(neg_zero, -nl::min(), neg_zero, false, 0);
+    TEST_SPECIAL(neg_zero, nl::max(), neg_zero, false, 0);
+    TEST_SPECIAL(neg_zero, -nl::max(), neg_zero, false, 0);
+
+    // fmod (+inf, y) == nl::quiet_NaN() plus invalid exception.
+    TEST_SPECIAL(inf, 3.0, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(inf, -1.1L, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(inf, 0.0, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(inf, neg_zero, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(inf, nl::denorm_min(), nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(inf, nl::min(), nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(inf, nl::max(), nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(inf, inf, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(inf, neg_inf, nl::quiet_NaN(), true, FE_INVALID);
+
+    // fmod (-inf, y) == nl::quiet_NaN() plus invalid exception.
+    TEST_SPECIAL(neg_inf, 3.0, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(neg_inf, -1.1L, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(neg_inf, 0.0, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(neg_inf, neg_zero, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(neg_inf, nl::denorm_min(), nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(neg_inf, nl::min(), nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(neg_inf, nl::max(), nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(neg_inf, inf, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(neg_inf, neg_inf, nl::quiet_NaN(), true, FE_INVALID);
+
+    // fmod (x, +0) == nl::quiet_NaN() plus invalid exception.
+    TEST_SPECIAL(3.0, 0.0, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(-1.1L, 0.0, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(0.0, 0.0, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(neg_zero, 0.0, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(nl::denorm_min(), 0.0, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(nl::min(), 0.0, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(nl::max(), 0.0, nl::quiet_NaN(), true, FE_INVALID);
+
+    // fmod (x, -0) == nl::quiet_NaN() plus invalid exception.
+    TEST_SPECIAL(3.0, neg_zero, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(-1.1L, neg_zero, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(0.0, neg_zero, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(neg_zero, neg_zero, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(nl::denorm_min(), neg_zero, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(nl::min(), neg_zero, nl::quiet_NaN(), true, FE_INVALID);
+    TEST_SPECIAL(nl::max(), neg_zero, nl::quiet_NaN(), true, FE_INVALID);
+
+    // fmod (x, +inf) == x for x not infinite.
+    TEST_SPECIAL(0.0, inf, 0.0, false, 0);
+    TEST_SPECIAL(neg_zero, inf, neg_zero, false, 0);
+    TEST_SPECIAL(nl::denorm_min(), inf, nl::denorm_min(), false, 0);
+    TEST_SPECIAL(nl::min(), inf, nl::min(), false, 0);
+    TEST_SPECIAL(nl::max(), inf, nl::max(), false, 0);
+    TEST_SPECIAL(3.0, inf, 3.0, false, 0);
+    // fmod (x, -inf) == x for x not infinite.
+    TEST_SPECIAL(0.0, neg_inf, 0.0, false, 0);
+    TEST_SPECIAL(neg_zero, neg_inf, neg_zero, false, 0);
+    TEST_SPECIAL(nl::denorm_min(), neg_inf, nl::denorm_min(), false, 0);
+    TEST_SPECIAL(nl::min(), neg_inf, nl::min(), false, 0);
+    TEST_SPECIAL(nl::max(), neg_inf, nl::max(), false, 0);
+    TEST_SPECIAL(3.0, neg_inf, 3.0, false, 0);
+
+    TEST_SPECIAL(0.0, nl::quiet_NaN(), nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(0.0, -nl::quiet_NaN(), nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(neg_zero, nl::quiet_NaN(), nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(neg_zero, -nl::quiet_NaN(), nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(1.0, nl::quiet_NaN(), nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(1.0, -nl::quiet_NaN(), nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(inf, nl::quiet_NaN(), nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(inf, -nl::quiet_NaN(), nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(neg_inf, nl::quiet_NaN(), nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(neg_inf, -nl::quiet_NaN(), nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(0.0, nl::signaling_NaN(), nl::quiet_NaN(), false, FE_INVALID);
+    TEST_SPECIAL(0.0, -nl::signaling_NaN(), nl::quiet_NaN(), false, FE_INVALID);
+    TEST_SPECIAL(neg_zero, nl::signaling_NaN(), nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(neg_zero, -nl::signaling_NaN(), nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(1.0, nl::signaling_NaN(), nl::quiet_NaN(), false, FE_INVALID);
+    TEST_SPECIAL(1.0, -nl::signaling_NaN(), nl::quiet_NaN(), false, FE_INVALID);
+    TEST_SPECIAL(inf, nl::signaling_NaN(), nl::quiet_NaN(), false, FE_INVALID);
+    TEST_SPECIAL(inf, -nl::signaling_NaN(), nl::quiet_NaN(), false, FE_INVALID);
+    TEST_SPECIAL(neg_inf, nl::signaling_NaN(), nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(neg_inf, -nl::signaling_NaN(), nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(nl::quiet_NaN(), 0.0, nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(-nl::quiet_NaN(), 0.0, nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(nl::quiet_NaN(), neg_zero, nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(-nl::quiet_NaN(), neg_zero, nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(nl::quiet_NaN(), 1.0, nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(-nl::quiet_NaN(), 1.0, nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(nl::quiet_NaN(), inf, nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(-nl::quiet_NaN(), inf, nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(nl::quiet_NaN(), neg_inf, nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(-nl::quiet_NaN(), neg_inf, nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(nl::signaling_NaN(), 0.0, nl::quiet_NaN(), false, FE_INVALID);
+    TEST_SPECIAL(-nl::signaling_NaN(), 0.0, nl::quiet_NaN(), false, FE_INVALID);
+    TEST_SPECIAL(nl::signaling_NaN(), neg_zero, nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(-nl::signaling_NaN(), neg_zero, nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(nl::signaling_NaN(), 1.0, nl::quiet_NaN(), false, FE_INVALID);
+    TEST_SPECIAL(-nl::signaling_NaN(), 1.0, nl::quiet_NaN(), false, FE_INVALID);
+    TEST_SPECIAL(nl::signaling_NaN(), inf, nl::quiet_NaN(), false, FE_INVALID);
+    TEST_SPECIAL(-nl::signaling_NaN(), inf, nl::quiet_NaN(), false, FE_INVALID);
+    TEST_SPECIAL(nl::signaling_NaN(), neg_inf, nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(-nl::signaling_NaN(), neg_inf, nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(nl::quiet_NaN(), nl::quiet_NaN(), nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(nl::quiet_NaN(), -nl::quiet_NaN(), nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(-nl::quiet_NaN(), nl::quiet_NaN(), nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(-nl::quiet_NaN(), -nl::quiet_NaN(), nl::quiet_NaN(), false, 0);
+    TEST_SPECIAL(nl::quiet_NaN(), nl::signaling_NaN(), nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(nl::quiet_NaN(), -nl::signaling_NaN(), nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(-nl::quiet_NaN(), nl::signaling_NaN(), nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(-nl::quiet_NaN(), -nl::signaling_NaN(), nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(nl::signaling_NaN(), nl::quiet_NaN(), nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(nl::signaling_NaN(), -nl::quiet_NaN(), nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(-nl::signaling_NaN(), nl::quiet_NaN(), nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(-nl::signaling_NaN(), -nl::quiet_NaN(), nl::quiet_NaN(), false,
+                 FE_INVALID);
+    TEST_SPECIAL(nl::signaling_NaN(), nl::signaling_NaN(), nl::quiet_NaN(),
+                 false, FE_INVALID);
+    TEST_SPECIAL(nl::signaling_NaN(), -nl::signaling_NaN(), nl::quiet_NaN(),
+                 false, FE_INVALID);
+    TEST_SPECIAL(-nl::signaling_NaN(), nl::signaling_NaN(), nl::quiet_NaN(),
+                 false, FE_INVALID);
+    TEST_SPECIAL(-nl::signaling_NaN(), -nl::signaling_NaN(), nl::quiet_NaN(),
+                 false, FE_INVALID);
+
+    TEST_SPECIAL(6.5, 2.25L, 2.0L, false, 0);
+    TEST_SPECIAL(-6.5, 2.25L, -2.0L, false, 0);
+    TEST_SPECIAL(6.5, -2.25L, 2.0L, false, 0);
+    TEST_SPECIAL(-6.5, -2.25L, -2.0L, false, 0);
+
+    TEST_SPECIAL(nl::max(), nl::max(), 0.0, false, 0);
+    TEST_SPECIAL(nl::max(), -nl::max(), 0.0, false, 0);
+    TEST_SPECIAL(nl::max(), nl::min(), 0.0, false, 0);
+    TEST_SPECIAL(nl::max(), -nl::min(), 0.0, false, 0);
+    TEST_SPECIAL(nl::max(), nl::denorm_min(), 0.0, false, 0);
+    TEST_SPECIAL(nl::max(), -nl::denorm_min(), 0.0, false, 0);
+    TEST_SPECIAL(-nl::max(), nl::max(), neg_zero, false, 0);
+    TEST_SPECIAL(-nl::max(), -nl::max(), neg_zero, false, 0);
+    TEST_SPECIAL(-nl::max(), nl::min(), neg_zero, false, 0);
+    TEST_SPECIAL(-nl::max(), -nl::min(), neg_zero, false, 0);
+    TEST_SPECIAL(-nl::max(), nl::denorm_min(), neg_zero, false, 0);
+    TEST_SPECIAL(-nl::max(), -nl::denorm_min(), neg_zero, false, 0);
+
+    TEST_SPECIAL(nl::min(), nl::max(), nl::min(), false, 0);
+    TEST_SPECIAL(nl::min(), -nl::max(), nl::min(), false, 0);
+    TEST_SPECIAL(nl::min(), nl::min(), 0.0, false, 0);
+    TEST_SPECIAL(nl::min(), -nl::min(), 0.0, false, 0);
+    TEST_SPECIAL(nl::min(), nl::denorm_min(), 0.0, false, 0);
+    TEST_SPECIAL(nl::min(), -nl::denorm_min(), 0.0, false, 0);
+    TEST_SPECIAL(-nl::min(), nl::max(), -nl::min(), false, 0);
+    TEST_SPECIAL(-nl::min(), -nl::max(), -nl::min(), false, 0);
+    TEST_SPECIAL(-nl::min(), nl::min(), neg_zero, false, 0);
+    TEST_SPECIAL(-nl::min(), -nl::min(), neg_zero, false, 0);
+    TEST_SPECIAL(-nl::min(), nl::denorm_min(), neg_zero, false, 0);
+    TEST_SPECIAL(-nl::min(), -nl::denorm_min(), neg_zero, false, 0);
+
+    TEST_SPECIAL(nl::denorm_min(), nl::max(), nl::denorm_min(), false, 0);
+    TEST_SPECIAL(nl::denorm_min(), -nl::max(), nl::denorm_min(), false, 0);
+    TEST_SPECIAL(nl::denorm_min(), nl::min(), nl::denorm_min(), false, 0);
+    TEST_SPECIAL(nl::denorm_min(), -nl::min(), nl::denorm_min(), false, 0);
+    TEST_SPECIAL(nl::denorm_min(), nl::denorm_min(), 0.0, false, 0);
+    TEST_SPECIAL(nl::denorm_min(), -nl::denorm_min(), 0.0, false, 0);
+    TEST_SPECIAL(-nl::denorm_min(), nl::max(), -nl::denorm_min(), false, 0);
+    TEST_SPECIAL(-nl::denorm_min(), -nl::max(), -nl::denorm_min(), false, 0);
+    TEST_SPECIAL(-nl::denorm_min(), nl::min(), -nl::denorm_min(), false, 0);
+    TEST_SPECIAL(-nl::denorm_min(), -nl::min(), -nl::denorm_min(), false, 0);
+    TEST_SPECIAL(-nl::denorm_min(), nl::denorm_min(), neg_zero, false, 0);
+    TEST_SPECIAL(-nl::denorm_min(), -nl::denorm_min(), neg_zero, false, 0);
+  }
+
+  void testRegularExtreme(FModFunc f) {
+
+    TEST_REGULAR(0x1p127L, 0x3p-149L, 0x1p-149L);
+    TEST_REGULAR(0x1p127L, -0x3p-149L, 0x1p-149L);
+    TEST_REGULAR(0x1p127L, 0x3p-148L, 0x1p-147L);
+    TEST_REGULAR(0x1p127L, -0x3p-148L, 0x1p-147L);
+    TEST_REGULAR(0x1p127L, 0x3p-126L, 0x1p-125L);
+    TEST_REGULAR(0x1p127L, -0x3p-126L, 0x1p-125L);
+    TEST_REGULAR(-0x1p127L, 0x3p-149L, -0x1p-149L);
+    TEST_REGULAR(-0x1p127L, -0x3p-149L, -0x1p-149L);
+    TEST_REGULAR(-0x1p127L, 0x3p-148L, -0x1p-147L);
+    TEST_REGULAR(-0x1p127L, -0x3p-148L, -0x1p-147L);
+    TEST_REGULAR(-0x1p127L, 0x3p-126L, -0x1p-125L);
+    TEST_REGULAR(-0x1p127L, -0x3p-126L, -0x1p-125L);
+
+    if constexpr (sizeof(T) >= sizeof(double)) {
+      TEST_REGULAR(0x1p1023L, 0x3p-1074L, 0x1p-1073L);
+      TEST_REGULAR(0x1p1023L, -0x3p-1074L, 0x1p-1073L);
+      TEST_REGULAR(0x1p1023L, 0x3p-1073L, 0x1p-1073L);
+      TEST_REGULAR(0x1p1023L, -0x3p-1073L, 0x1p-1073L);
+      TEST_REGULAR(0x1p1023L, 0x3p-1022L, 0x1p-1021L);
+      TEST_REGULAR(0x1p1023L, -0x3p-1022L, 0x1p-1021L);
+      TEST_REGULAR(-0x1p1023L, 0x3p-1074L, -0x1p-1073L);
+      TEST_REGULAR(-0x1p1023L, -0x3p-1074L, -0x1p-1073L);
+      TEST_REGULAR(-0x1p1023L, 0x3p-1073L, -0x1p-1073L);
+      TEST_REGULAR(-0x1p1023L, -0x3p-1073L, -0x1p-1073L);
+      TEST_REGULAR(-0x1p1023L, 0x3p-1022L, -0x1p-1021L);
+      TEST_REGULAR(-0x1p1023L, -0x3p-1022L, -0x1p-1021L);
+    }
+  }
+};
+
+#define LIST_FMOD_TESTS(T, func)                                               \
+  using LlvmLibcFmodTest = FmodTest<T>;                                        \
+  TEST_F(LlvmLibcFmodTest, SpecialNumbers) { testSpecialNumbers(&func); }      \
+  TEST_F(LlvmLibcFmodTest, RegularExtreme) { testRegularExtreme(&func); }
+
+#endif // LLVM_LIBC_TEST_SRC_MATH_FMODTEST_H

libc/test/src/math/differential_testing/CMakeLists.txt

add_diff_binary(
   SRCS
     hypot_perf.cpp
   DEPENDS
     .binary_op_single_output_diff
     libc.src.math.hypot
   COMPILE_OPTIONS
     -fno-builtin
 )
+
+add_diff_binary(
+  fmodf_diff
+  SRCS
+    fmodf_diff.cpp
+  DEPENDS
+    .single_input_single_output_diff
+    libc.src.math.fmodf
+)
+
+add_diff_binary(
+  fmodf_perf
+  SRCS
+    fmodf_perf.cpp
+  DEPENDS
+    .single_input_single_output_diff
+    libc.src.math.fmodf
+  COMPILE_OPTIONS
+    -fno-builtin
+)
+
+add_diff_binary(
+  fmod_diff
+  SRCS
+    fmod_diff.cpp
+  DEPENDS
+    .single_input_single_output_diff
+    libc.src.math.fmod
+)
+
+add_diff_binary(
+  fmod_perf
+  SRCS
+    fmod_perf.cpp
+  DEPENDS
+    .single_input_single_output_diff
+    libc.src.math.fmod
+  COMPILE_OPTIONS
+    -fno-builtin
+)

libc/test/src/math/differential_testing/fmod_diff.cpp

+//===-- Differential test for fmod ----------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "BinaryOpSingleOutputDiff.h"
+
+#include "src/math/fmod.h"
+
+#include <math.h>
+
+BINARY_OP_SINGLE_OUTPUT_DIFF(double, __llvm_libc::fmod, ::fmod, "fmod_diff.log")

libc/test/src/math/differential_testing/fmod_perf.cpp

+//===-- Differential test for fmod ----------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "BinaryOpSingleOutputDiff.h"
+
+#include "src/math/fmod.h"
+
+#include <math.h>
+
+BINARY_OP_SINGLE_OUTPUT_PERF(double, __llvm_libc::fmod, ::fmod, "fmod_perf.log")

libc/test/src/math/differential_testing/fmodf_diff.cpp

+//===-- Differential test for fmodf ---------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "BinaryOpSingleOutputDiff.h"
+
+#include "src/math/fmodf.h"
+
+#include <math.h>
+
+BINARY_OP_SINGLE_OUTPUT_DIFF(float, __llvm_libc::fmodf, ::fmodf,
+                             "fmodf_diff.log")

libc/test/src/math/differential_testing/fmodf_perf.cpp

+//===-- Differential test for fmodf ---------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "BinaryOpSingleOutputDiff.h"
+
+#include "src/math/fmodf.h"
+
+#include <math.h>
+
+BINARY_OP_SINGLE_OUTPUT_PERF(float, __llvm_libc::fmodf, ::fmodf,
+                             "fmodf_perf.log")

libc/test/src/math/exhaustive/CMakeLists.txt

add_fp_unittest(
   DEPENDS
     .exhaustive_test
     libc.include.math
     libc.src.math.hypotf
     libc.src.__support.FPUtil.fputil
   LINK_LIBRARIES
     -lpthread
 )
+
+add_fp_unittest(
+  fmod_generic_impl_test
+  NO_RUN_POSTBUILD
+  NEED_MPFR
+  SUITE
+    libc_math_exhaustive_tests
+  SRCS
+    fmod_generic_impl_test.cpp
+  DEPENDS
+    libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.generic.fmod
+)

libc/test/src/math/exhaustive/fmod_generic_impl_test.cpp

+//===-- Utility class to test FMod generic implementation -------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+#include "src/__support/CPP/TypeTraits.h"
+#include "src/__support/FPUtil/generic/FMod.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/FPMatcher.h"
+#include "utils/UnitTest/Test.h"
+
+#include <array>
+#include <limits>
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+template <typename T, bool InverseMultiplication>
+class LlvmLibcFModTest : public __llvm_libc::testing::Test {
+
+  using DivisionHelper = __llvm_libc::cpp::ConditionalType<
+      InverseMultiplication,
+      __llvm_libc::fputil::generic::FModDivisionInvMultHelper<T>,
+      __llvm_libc::fputil::generic::FModDivisionSimpleHelper<T>>;
+
+  static constexpr std::array<T, 11> test_bases = {
+      T(0.0),
+      T(1.0),
+      T(3.0),
+      T(27.0),
+      T(11.0 / 8.0),
+      T(2.764443),
+      T(1.0) - std::numeric_limits<T>::epsilon(),
+      T(1.0) + std::numeric_limits<T>::epsilon(),
+      T(M_PI),
+      T(M_SQRT2),
+      T(M_E)};
+
+public:
+  void testExtensive() {
+    using FMod = __llvm_libc::fputil::generic::FMod<
+        T, __llvm_libc::fputil::generic::FModFastMathWrapper<T>,
+        DivisionHelper>;
+    using nl = std::numeric_limits<T>;
+    int min2 = nl::min_exponent - nl::digits - 5;
+    int max2 = nl::max_exponent + 3;
+    for (T by : test_bases) {
+      for (int iy = min2; iy < max2; iy++) {
+        T y = by * std::ldexp(2, iy);
+        if (y == 0 || !std::isfinite(y))
+          continue;
+        for (T bx : test_bases) {
+          for (int ix = min2; ix < max2; ix++) {
+            T x = bx * std::ldexp(2, ix);
+            if (!std::isfinite(x))
+              continue;
+            T result = FMod::eval(x, y);
+            mpfr::BinaryInput<T> input{x, y};
+            EXPECT_MPFR_MATCH(mpfr::Operation::Fmod, input, result, 0.0);
+          }
+        }
+      }
+    }
+  }
+};
+
+using LlvmLibcFModFloatTest = LlvmLibcFModTest<float, false>;
+TEST_F(LlvmLibcFModFloatTest, ExtensiveTest) { testExtensive(); }
+
+using LlvmLibcFModFloatInvTest = LlvmLibcFModTest<float, true>;
+TEST_F(LlvmLibcFModFloatInvTest, ExtensiveTest) { testExtensive(); }
+
+using LlvmLibcFModDoubleTest = LlvmLibcFModTest<double, false>;
+TEST_F(LlvmLibcFModDoubleTest, ExtensiveTest) { testExtensive(); }
+
+using LlvmLibcFModDoubleInvTest = LlvmLibcFModTest<double, true>;
+TEST_F(LlvmLibcFModDoubleInvTest, ExtensiveTest) { testExtensive(); }

libc/test/src/math/fmod_test.cpp

+//===-- Unittests for fmod ------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "FModTest.h"
+
+#include "src/math/fmod.h"
+
+LIST_FMOD_TESTS(double, __llvm_libc::fmod)

libc/test/src/math/fmodf_test.cpp

+//===-- Unittests for fmodf -----------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "FModTest.h"
+
+#include "src/math/fmodf.h"
+
+LIST_FMOD_TESTS(float, __llvm_libc::fmodf)

libc/utils/MPFRWrapper/MPFRUtils.h

enum class Operation : int {
   BeginUnaryOperationsTwoOutputs,
   Frexp, // Floating point output, the first output, is the fractional part.
   EndUnaryOperationsTwoOutputs,
 
   // Operations wich take two floating point nubmers of the same type as
   // input and produce a single floating point number of the same type as
   // output.
   BeginBinaryOperationsSingleOutput,
+  Fmod,
   Hypot,
   EndBinaryOperationsSingleOutput,
 
   // Operations which take two floating point numbers of the same type as
   // input and produce two outputs. The first output is a floating nubmer of
   // the same type as the inputs. The second output is af type 'int'.
   BeginBinaryOperationsTwoOutputs,
   RemQuo, // The first output, the floating point output, is the remainder.

libc/utils/MPFRWrapper/MPFRUtils.cpp

public:
   }
 
   MPFRNumber floor() const {
     MPFRNumber result(*this);
     mpfr_floor(result.value, value);
     return result;
   }
 
+  MPFRNumber fmod(const MPFRNumber &b) {
+    MPFRNumber result(*this);
+    mpfr_fmod(result.value, value, b.value, mpfr_rounding);
+    return result;
+  }
+
   MPFRNumber frexp(int &exp) {
     MPFRNumber result(*this);
     mpfr_exp_t resultExp;
     mpfr_frexp(&resultExp, result.value, value, mpfr_rounding);
     exp = resultExp;
     return result;
   }
 

 
 template <typename InputType>
 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
 binary_operation_one_output(Operation op, InputType x, InputType y,
                             unsigned int precision, RoundingMode rounding) {
   MPFRNumber inputX(x, precision, rounding);
   MPFRNumber inputY(y, precision, rounding);
   switch (op) {
+  case Operation::Fmod:
+    return inputX.fmod(inputY);
   case Operation::Hypot:
     return inputX.hypot(inputY);
   default:
     __builtin_unreachable();
   }
 }
 
 template <typename InputType>